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The  primary  purpose  of  this  research  was  to  investigate  possible 
means  of  improving  imaging  of  breast  cancers  by  transillumination  using 
near  infrared  radiation.  The  need  for  such  a  detection  scheme  was  clear 
considering  the  possible  carcenogenic  properties  of  x-ray  mammography 
which  is  the  standard  modality  for  breast  cancer  detection  today. 

Lack  of  experimental  data  related  to  optical  parameters  of  breast 
tissue  set  limits  on  the  degree  of  applicability  of  the  methods 
presented.  Nevertheless,  the  model  of  light  propagation  in  the  breast 
employed  in  the  study  yielded  results  for  transmittance  which  are  in 
agreement  with  actual  breast  trans illumination  to  a  remarkable  degree 
considering  the  paucity  of  quantitative  data  available. 

The  encouraging  results  obtained  with  time  gating  and  spatial 
filtering  (described  in  Chapter  III)  of  light  propagating  through  a 
dilute  blood  medium  should  be  followed  up  by  applying  the  same  method  to 
breast  trans illumination. 
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Abstract 

This  investigation  is  a  theoretical  analysis  of  the  propagation  of 
near  infrared  light  through  human  breast  tissue  with  a  view  to 
improving  imaging  of  breast  cancers  by  transillumination.  The  analysis 
employed  multiple  scattering  and  radiation  transport  theory. 

Calculation  showed  that  breast  tissue  is  optically  thick  and  does 
not  allow  for  imaging  details  of  tumors  which  may  be  imbedded  in  the 
tissue.  An  example  calculation  carried  out  on  the  method  of  phase 
retrieval  demonstrated  that  the  error  in  the  estimate  of  the  Fourier 
modulus  is,  for  normal  breast  tissue,  nearly  100  percent  and  that, 
therefore,  the  phase  of  an  object  (tumor)  cannot  be  retrieved  by  this 
method. 

A  Monte  Carlo  simulation  employing  time  gating  and  spatial  filtering 
(TGSF)  was  presented  showing  that,  in  some  cases,  an  improvement  in 
image  contrast  of  77  percent  is  theoretically  possible  when 
transilluminating  a  dilute  blood  medium  containing  a  blood  vessel.  A 
calculation  of  total  transmittance  obtained  by  Monte  Carlo  simulation  on 
one  hand  and  by  multiple  scattering  and  radiation  transport  theory  on 
the  other,  indicated  an  essential  agreement  of  the  results  for  tissue 
thicknesses  less  than  60mm. 

A  method  of  tying  the  results  of  Monte  Carlo  simulation  to  those  of 
multiple  scattering  and  radiation  transport  is  presented.  An  effective 
optical  thickness  reduced  by  TGSF  is  defined  and  a  method  is  described 
to  determine  if  the  reduced  effective  thickness  allows  for  reasonable 
images  to  be  obtained. 
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THEORETICAL  ANALYSIS  OF  CANCER  DETECTION 
IN  THE  HUMAN  BREAST  BY  TRANSILLUMINATION 


I.  Trans illumination 


Introduction 

Concern  over  the  limitations  and  risks  (6)  associated  with 
x-radiation  in  the  imaging  and  diagnosis  of  human  breast  cancer  is 
motivating  researchers  to  seek  other,  safer  methods  of  detecting  breast 
lesions.  One  such  method,  trans illumination,  is  a  technique  of  imaging 
structures  of  human  tissue  by  propagating  light  through  them  (See  Figure 
1).  The  technique  has  gone  by  different  names  depending  on  light  sources 
and  methods  of  light  detection  used.  The  term  transillumination 
initially  served  to  name  a  procedure  which  used  white  light  as  a  source 
and  where  observation  of  the  illuminated  tissue  was  made  by  the  unaided 
eye.  In  cases  where  white  light  sources  are  used,  but  detection  is 
performed  via  infrared  photography,  the  trans illumination  procedure  has 
been  called  diaphanography  (23:123).  The  term  diaphanography  will  not, 
however,  be  used  in  this  paper  as  the  more  advanced  transillumination 
methods  employ  infrared  illumination  (by  filtering  white  light)  and 
infrared  video  camera  detection,  amplification,  and  video  display 
schemes.  In  these  cases,  transillumination  is  referred  to  as  light 
scanning  or  infrared  light  scanning  (5:409).  The  reader  will  note  that 
the  term  transillumination  can  be  used  at  times,  in  its  generic  sense, 
to  refer  to  any  of  the  optical  imaging  techniques  mentioned  above. 

Structures  within  the  human  breast  have  varying  optical  properties 


such  that,  when  the  structures  are  trans illuminated,  variations  result 
in  the  distribution  of  light  intensity  and  color  in  their  images. 
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Figure  1.  (a):  Carcinoma  imaged  by  transillumination  (23:127), 

(b):  Same  carcinoma  imaged  by  x-ray  manmography  (23:127). 


Researchers  have  found  that  breast  cancers  may  be  identified  by  the 
shadows  they  cast,  and  benign  lesions  by  the  areas  of  increased 
brightness  they  create  in  the  image  (39:14).  However,  the 
transillunination  technique  has,  to  date,  demonstrated  only  sufficient 
effectiveness  to  perform  as  an  augmentation  to  traditional  marnnography 
and  has  not  reached  a  level  of  development  whereby  it  could  be 
considered  as  an  alternative  to  mammography  as  the  primary  imaging  tool 
in  the  diagnosis  of  breast  cancer. 

The  following  section  of  Chapter  I  is  a  review  of  the  current 
literature  relevant  to  transillumination  as  an  aid  in  the  diagnosis  of 
diseased  breast  tissue.  The  references  cited  were  collected  by  means  of 
a  computerized  literature  search  using  the  medical  library  of  DIALOG 
software.  The  reader  will  find  that,  save  for  the  article  by  Watmough 
(which  is  the  last  presented  here),  the  review  is  organized 
chronologically,  showing  the  improvements  made,  and  as  a  consequence, 
the  problems  which  have  confronted  and  continue  to  define  the  limits  of 
the  applicability  of  transillunination. 

Background 

The  first  diagnostic  use  of  transillumination  of  the  breast  was  made 
by  Culter  in  1929  (9).  He  employed  visible  wavelengths  of  white  light 
and  reported  some  success  in  identifying  tissue  variations.  However,  it 
was  recognized  early  on  that  "it  was  not  possible  to  differentiate 
shadows  of  benign  from  malignant  lesions,  and  that  light  diffusion  from 
the  available  probe  often  prevented  smaller  lesions  from  being  seen" 
(23:123). 


Lafreniere,  et  al.  reported  that,  in  spite  of  initial  support  for 
the  new  transillumination  technique,  the  procedure  fell  into  disrepute, 
probably  due  to  inadequacy  of  the  available  equipment  (23:123).  Over 


the  years  a  number  of  improvements  were  made  to  both  the  equipment  and 
technique.  Among  these  improvements  were  (1)  illumination  with  longer 
wavelength,  narrower  spectrum  sources  which  increased  tissue 
penetration  and  reduced  scattering,  (2)  the  introduction  of  infrared 
photography,  and  (3)  the  use  of  fiber  optic  bundles  to  create  a  cold 
light  source  to  convey  the  light  energy  to  the  tissue  (this  overcame  an 
early  problem  with  heat  generated  by  the  light  source).  Despite  these 
improvements,  reports  indicated  the  technique  was  cumbersome  and  did  not 
add  significantly  to  the  conventional  physical  examination  and  mammogram 
(18). 

Better  results  were  obtained  by  Morton  (25)  who,  in  1981,  reported 
using  a  red  and  infrared  light  source  to  illuminate  the  tissue.  In 
addition,  he  introduced  an  infrared  video  camera  and  displayed  its 
images  on  a  black  and  white  video  monitor.  This  allowed  for 
amplification  of  the  image  and  thereby  reduced  the  intensity  requirement 
of  the  source.  As  a  result  of  Morton's  improvements,  internal 
structures  which  had  previously  been  unrecognizable  to  the  unaided  eye 
could  be  viewed  and  identified  on  the  video  monitor  (26). 

Morton's  technique  had  been  based  on  Cartwright's  (7)  observation 
made  in  1930  that  greater  transmission  through  human  tissue  was  possible 
with  infrared  illumination.  The  work  of  Anderson  and  Parrish  (2:17) 
supports  Cartwright's  findings  showing  that  transmittance  of  light 
through  human  skin  is  decreased  in  the  shorter  wavelength  spectrum  and 
increases  in  the  near  infrared. 

In  1983,  Bartrum  and  Crow  (5:409)  reported  that  transillumination 
light  scanning  can  yield  clinical  results  comparable  to  mammography  if 
the  personnel  who  are  interpreting  light  scan  images  are  trained  and 
have  some  knowledge  of  other  clinical  data  related  to  the  patient.  When 


the  light  scan  interpreter  is  blind  to  other  relevant  clinical  data,  a 
diagnosis  based  on  the  light  scan  is  generally  poorer  than  that  obtained 
from  mammography.  In  their  report,  Bartrum  and  Crow  cite  conditions 
they  feel  must  be  met  to  obtain  optimal  results  from  light  scanning. 
First,  diffuse  incident  light  is  required  since  breast  tissue  is  a 
highly  scattering  medium.  When  focused  light  is  used  the  light  scan 
becomes  more  difficult  to  interpret.  Second,  the  optical  absorption 
differences  between  cancerous  and  adjacent  tissue  can  be  blurred  if  the 
transilluminating  light  is  too  intense.  Third,  vigorous  physical 
compression  of  the  breast  is  necessary  since  the  farther  a  lesion  is 
located  from  the  skin  surface,  the  weaker  is  the  shadow  cast  by  the 
lesion.  Bartrum  and  Crow  found  that  the  "shadow  from  a  1-cm-diameter 
lesion  (in  a  breast  phantom)  begins  to  fade  rapidly  when  the  lesion  is 
more  than  15mm  from  the  skin  surface;  at  25mm,  the  shadow  is  obscured" 
(5:413).  For  this  reason  a  breast  lesion  may  be  invisible  if  "it  is  not 
compressed  to  within  2cm  of  a  skin  surface"  (5:413).  Fourth,  because  a 
lesion  may  become  undetectable  a  short  distance  beneath  the  tissue 
surface,  light  scanning  from  multiple  angles  is  required  to  ensure 
maximum  possible  coverage  of  the  tissue.  Fifth,  and  lastly,  since 
carcinomas  absorb  light  differentially  in  narrow  spectra,  a  narrow 
spectrum  of  illuminating  light  should  be  used  as  a  source  (this  could  be 
accomplished  by  filtering  a  white  light  source).  Bartrum  and  Crow  go  on 
to  say  ".  .  .we  performed  in  vivo  spectral  analysis  of  11  carcinomas.  .  . 
which  showed  preferential  tumor  absorption  in  two  relatively  narrow 
wavelenght  bands,  one  in  the  far  visible  red  spectrun  and  another  in  the 
infrared  spectrum"  (5:413).  In  1986,  Lafreniere  et  al.  (23)  repjorted 
the  results  of  a  clinical  evaluation  of  transillumination  and  infrared 
light  scanning  as  compered  with  conventional  physical  examination  and 


mammography.  Ihe  test  involved  600  patients  with  79  biopsies  performed 

revealing  26  carcinomas  (23:124).  These  researchers  reported  that 

examination,  conventional  manmography;  simple  transillumination, 
I.L.S.  (infrared  light  scanning),  and  simple  transillumination 
combined  with  I.L.S.  respectively  showed  a  true  positive 
interpretation  in  21,  22,  19,  25,  and  25  of  26  breast  lesions 
histologically  positive  for  malignancy  (sensitivity  81%,  85%,  73%, 
96%,  96%)  (19:124). 

These  results  indicate  that  infrared  light  scanning  missed  only  one 
of  26  carcinomas  verified  histologically.  In  this  case,  the 
interpreters  of  the  images  of  manmography,  transillumination,  or 
infrared  light  scanning  techniques  were  trained  in  each  technique  but 
had  no  knowledge  of  other  clinical  data  relating  to  the  patients. 

An  explanation  of  the  physical  mechanism  responsible  for  the 
appearance  of  the  images  produced  by  transillumination  and  infrared 
light  scanning  was  given  by  Watmough  in  1982  (39:14).  Watmough  had  set 
out  to  measure  the  transmission  coefficient  of  sections  of  tissue 
removed  from  the  female  breast  during  surgical  operation.  He  had 
expected  to  find  that  "different  lesions  would  transmit  light  to  the 
extent  dependent  on  the  typo  of  disease  process"  (39:11).  Instead,  he 
found  that  despite  large  variations  (size  and  shapo)  in  tissue  samples 
taken,  there  was  no  clear  distinction  between  transmission 
characteristics  of  malignant  tumors  and  benign  lesions.  Watmough  argued 
that  since  the  variations  in  light  transmission  were  not  caused  by  the 
diseased  tissues  themselves,  transillumination  and  light  scanning 
were  simply  not  imaging  the  diseased  tissues.  He  reports,  ".  .  .the 
spectro-photometric  traces  in  every  case  demonstrated  the  absorption 
bands  of  oxyhaemoglobin  [oxygenated  hemoglobin]  .  .  ."  (39:11).  From 
these  findings,  Watmough  postulated  that  the  transillumination  procedure 
actually  measured  the  number  density  of  red  cells  in  the  illuminated 
tissue.  To  determine  the  accuracy  of  his  hypothesis,  Watmough 
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trans illuminated  (with  a  xenon  flash  tube)  samples  of  whole  blood  at 
varying  degrees  of  dilution  and  exposed  infrared  color  film  to  the  light 
scan  images  of  the  breast.  A  spectral  analysis  of  samples  of  diluted 
blood  was  then  carried  out;  the  results  showing  that  the  samples  ceased 
to  transmit  light  at  frequencies  (colors)  related  to  the  color  of  whole 
blood  in  a  sample.  For  samples  with  higher  blood  densities  the 
frequency  cut-offs  occurred  in  the  longer  wavelength  regions  (red  and 
orange).  The  more  diluted  samples  showed  frequency  cut-offs  toward  the 
blue  end  of  the  spectrum.  Watmough  writes,  ".  .  .These  experiments 
provide  excellent  confirmation  of  how  red  cell  densities  in  effect 
provide  a  range  of  colour  filters.  .  .  ."(39:13). 

Breast  Anatomy 


The  Normal  Breast. 

The  description  given  here  of  the  anatomy  and  histology  of  the  human 
breast  is  taken  from  Egan  (12:17-24).  The  structure  of  the  normal 
breast  will  be  presented  first  to  be  followed  by  the  formulation  of  an 
optical  model  of  the  breast. 

The  breast  is  located  between  the  superficial  and  deep  layers  of  the 
superficial  fascia  (see  Figure  2).  The  superficial  layer  of  fascia, 
seperated  from  the  skin  by  0.5  to  2.5  cm  of  subcutaneous  fat  and  areolar 
tissue  (the  areola  is  the  region  of  darker  skin  surrounding  the  nipple), 
forms  an  irregular  boundary  for  the  surface  of  the  glandular  tissue. 
Projections  of  collagenous  fiberous  tissue  extend  from  the  fascia  and 
surround  lobules  of  fat.  These  masses  of  collagen  fibers  (low  in  water 
content)  form  suspensory  ligaments  called  Cooper’s  ligaments  (12:17). 
Younger  women  tend  to  have  breast  tissue  which  is  primarily  fibrous  with 
comparatively  small  amounts  of  fat  (28).  This  fact  is  of  particular 
importance  since  transillumination  is  targeted  at  younger  persons  below 
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Figure  2.  The  Normal  Breast  (12:18).  An  artist's  conception  of  the 
various  anatomic  structures  of  the  breast. 
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the  age  of  50  (see  Table  I).  This  means  that,  for  the  purposes  of 
transillunination,  the  breast  may  be  considered  primarily  fibrous.  When 
fat  lobules  occur  they  do  so  with  dimensions  on  the  order  of 
centimeters  (8:43-48). 

Within  the  breast  are  fifteen  to  twenty  irregular  lobes,  converging 
to  the  nipple.  Each  lobe  is  drained  by  its  own  lactiferous  duct,  and  is 
provided  with  an  excretory  duct,  2  to  4.5  run  in  diameter  (12:17).  These 
lobes  are  macroscopic  in  size  and  are,  therefore,  many  times  larger  than 
the  collagen  fibers  which  have  dimensions  on  the  order  of  a  hundred 
nanometers. 

Ihe  anatomic  structures  of  the  breast  are,  for  purposes  of 
classification,  divided  between  functional  structures  called  the 
parenchyma,  and  supporting  structures  called  the  stroma  (see  Figure  3). 
Because  the  parenchyma  (such  as  the  lactiferous  ducts  and  lobes)  tend  to 
be  very  large  compared  to  wavelength  of  radiation  in  transillumination, 
they  cannot  be  considered  to  contribute  greatly  to  light  scattering- 
also,  in  the  non-lac tating  breast  they  are  of  low  water  content  and, 
therefore,  are  not  highly  absorbing.  For  these  reasons,  the  description 
given  here  of  the  macroscopic  structure  of  the  breast  concentrates  on 
the  fibrous  tissue  of  the  stroma.  The  stroma  consists  of  dense, 
collagenous,  intralobular  connective  tissue  (see  Figure  3),  containing 
the  large  blood  vessels,  nerves,  lymphatics  and  varying  amounts  of 
adipose  tissue.  This  connective  tissue  produces  the  poorly  defined 
septa  between  the  lobules  and  lobes  of  the  mammary  gland  (12:19). 

In  the  craniocaudad  view  (see  Figure  4)  the  dense  inner  triangle  of 
tissue  produces  the  varying  degrees  of  diffuse  opasity  to  x-rays. 
Contrast  (in  x-ray  manmography )  in  the  image  is  dependent  upon  the 
relative  density  of  fatty  and  fibrous  tissues.  A  reticular  appearance 

-9- 


■  .v*  r  ’ 

«  4  '  »  *  '  4 


TABLE  I 


BREAST  TYPES 

Classification  of  breast  types  showing  percent  fibrous  and  fatty 
tissues  for  various  age  groups. 

FIBROUS  FATTY 


Age  Group 

Number 

Percent 

Number 

Percent 

Total 

21-25 

13 

100 

- 

- 

13 

26-30 

11 

91.6 

1 

8.4 

12 

31-35 

39 

88.7 

5 

11.3 

44 

36-40 

27 

81.8 

6 

11.2 

33 

41-45 

34 

75.5 

11 

24.5 

45 

46-50 

38 

65.5 

20 

34.5 

58 

51-55 

54 

55.1 

44 

44.9 

98 

56-60 

35 

37.2 

59 

62.8 

94 

61-65 

55 

39.8 

83 

60.2 

138 

66-70 

41 

35.9 

73 

64.1 

114 

71-75 

51 

52.2 

46 

47.4 

97 

76-80 

24 

64.8 

13 

35.2 

37 

81-85 

8 

66.6 

4 

33.4 

12 

86-90 

4 

80.0 

1 

20.0 

5 

21-90 

434 

54.3 

366 

45.8 

800 

(28 
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Figure  4.  Normal  Craniocaudad  View  of  the  Breast  (12:24).  A 
thirty-one  year  old  white  female  with  mammogram  done  for  routine 
check-up. 


caused  by  the  stroma  or  connective  tissue  framework  is  evident  in  x-ray 
images;  the  connective  tissue  may  be  nearly  homogeneous  in  the  very 
dense  breast  (12:22).  Again,  in  transillumination,  the  fatty  lobules 
are  not  considered  either  highly  scattering  or  absorbing  (31). 

To  the  extent  that  the  breast  is  a  light  absorbing  optical 
mediim,  the  absorption  experienced  is  due,  on  a  comparative  basis, 
almost  entirely  to  the  freely  flowing  plasma  in  veins  and  arteries. 
Accurate  values  for  the  volumetric  content  of  blood  in  the  breast  were 
not  available  for  this  writing.  The  subject  of  the  scattering  and 
absorption  characteristics  of  whole  blood  (blood  cells  and  plasma)  will 
be  covered  in  greater  detail  in  the  theoretical  model  of  breast  to  be 


described  below. 

An  Optical  Model  for  the  Human  Breast. 

Theoretical  models  of  the  breast  where  scattering  is  of  relative 
unimportance  often  consider  the  breast  to  be  a  medium  made  up  of  water 
and  fat  (27).  This  approximation  seems  suited  to  irradiation  of  the 
breast  with  x-rays  where  absorption  plays  the  primary  role  in 
attenuation.  However,  in  the  red  and  near  infrared  region  of  the 
spectrum  where  transillumination  occurs,  both  absorption  and  scattering 
play  a  role  in  determining  transmission  properties  of  the  breast. 

Absorption  Properties  of  the  Breast. 

By  far,  the  greatest  amount  of  absorption  in  biological  media  is 
found  to  occur  in  blood-borne  pigments  and  blood  plasma  (2:17).  In 
tissues  of  low  blood  or  water  content  relatively  little  absorption  is 
noted  (21:695).  Specific  optical  data  on  the  breast  which  would  allow 
one  to  isolate  transmission  losses  due  to  scattering  from  those  due  to 
absorption  are  lacking.  The  optical  characteristics  of  the  dermis  of 
the  skin  are,  however,  in  this  research  considered  to  have  an  especially 


close  affinity  to  the  optical  properties  of  the  breast  as  the  dermis  has 
a  high  content  of  collagenous  fibers  (these  fibers  account  for  70 
percent  of  the  dry  weight  of  the  dermis)  (3:256).  Anderson  and  Parrish 
state  (2:17)  with  regard  to  Findlay's  (16)  optical  measurements  of  pig 
dermis,  "Sunming  Findlay's  transmittance  and  remittance  spectra  gives 
values  close  to  1.0  (100  percent)  across  the  entire  visible  spectrum, 
indicating  that  very  little  visible  light  was  actually  absorbed." 

Figure  5  shows  spectral  transmittance  and  remittance  data  (analagous  to 
those  of  Findlay)  for  a  200  m  thick  section  of  human  dermis 
in  vitro.  Figure  6  presents  calculated  values  of  the  scattering 
coefficient  S  and  the  absorption  coefficient  K  for  human  dermis  in  vitro. 
The  coefficients  S  and  K  are  calculated  from  the  Kubelka-Munk  theory  of 
radiation  transfer  (22:106-116)  and  should  not  be  confused  with  the 
scattering  cross  section  (J  or  the  absorption  cross  section  (J 
which  will  be  used  in  Chapter  II  in  the  development  of  the  analytic  and 
transport  theories  of  scattering.  The  Kubelka-Munk  theory  is  not 
employed  in  this  research  as  it  does  not  lend  itself  readily  to 
interpretations  relating  to  actual  physical  parameters  of  individual 
particles.  This  should  not  be  construed,  however,  to  imply  that  the 
Kubelka-Munk  theory  does  not  give  valid  quantitative  results  (22:106). 

The  data  in  Figure  6  indicate  that  absorption  is  indeed  low  over 
the  visible  spectrum  and  also  in  the  near  infrared  region  except  at  the 
prominent  absorption  bands  of  water.  Figure  5  shows  peak  transmission 
occurs  at  about  1.23fJ,m  where,  again,  there  is  little  absorption. 

From  the  data  presented,  it  is  apparent  that  in  the  red  and  near 
infrared  portion  of  the  spectnm,  a  mediim  with  a  high  percentage  of 
collagenous  fibers  can  be  approximated  as  one  having  a  low  value  of 
absorptance.  Moreover,  due  to  the  fact  that  fatty  tissues  are  low  in 
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Figure  5  (2:17).  Spectral  Transmittance  and  Remittance  of  a  200  m 
Thick  Section  of  Human  Dennis.  Data  is  given  in  percent  remittance  (R) 
or  transmittance  (T). 
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Figure  6.  (2:17).  Diffuse  Scattering  (S)  and  Absorption  (K) 
Coefficients  for  Human  Dermis.  Data  is  calculated  from  the  Kubelka-Munk 
theory  of  radiation  transfer. 


water  content  (21:695),  one  may  conclude  that  these  tissues,  too, 
evidence  relatively  low  absorption.  It  is  not  possible  to  determine  the 


absolute  absorptance  of  the  fatty  tissues  at  near  Ijjm  wavelengths  (data 
in  the  literature  seems  to  be  concentrated  in  the  microwave  or 
ultraviolet  region  of  the  spectrum).  Calculations  to  be  done  in  Chapter 
III  will  show  that  total  theoretical  transmission  based  on  a  model 
assuming  low  absorptance  is  not  at  great  variance  with  transmission  data 
collected  under  actual  transillumination.  Consequently,  the  model  of 
breast  tissue  adopted  in  this  research  considers  breast  tissues  to 
exhibit  low  absorption  of  near  infrared  radiation  except  in  areas  where 
blood  plasma  is  present. 

Light  Scattering  Properties  of  the  Breast. 

In  this  section,  the  scattering  properties  of  fatty  tissues, 
collagen  fibers,  and  whole  blood  are  described. 

However,  before  considering  particular  tissue  types,  it  is 
convenient  to  describe  some  qualitative  aspects  of  scattering,  in 
general . 

First,  if  the  size  of  a  particle  is  much  smaller  than  a  wavelength, 

the  scattering  cross  section  (j  (a  measure  of  the  strength  of 

s 

scattering)  is  inversely  proportional  to  the  fourth  power  of  the 
wavelength  and  directly  proportional  to  the  square  of  the  volume  of  the 
particle  (19:13)  (see  Appendix  A  for  a  definition  ofCO.  This  type  of 
scattering  is  called  Rayleigh  scattering.  It  is  generally  isotropic 
(equal  in  all  directions)  in  its  radiation  pattern  and  comparatively 
weak  due  to  the  A”^  dependence.  As  particle  size  increases  and 
becomes  comparable  to  a  wavelength,  the  scattering  is  much  stronger 
than  in  Rayleigh  scattering  and  is  forward  directed  (or  in  a  dipole 
pattern).  When  particle  size  greatly  exceeds  the  wavelength,  the 


scattering  cross  section  is  proportional  to  X  ^  and  is  stronger  than 
Rayleigh  scattering  but  weaker  than  scattering  from  particles  close  to  a 
wavelength  in  size.  In  this  case,  the  scattering  is  more  forward 
directed.  This  third  classification  of  scattering  will  be  referred  to 
as  large  particle  scattering. 

Parrish,  in  discussing  the  types  of  scattering  (as  described  above) 
which  occur  in  skin  stated,  ".  .  .all  these  general  types  of  scattering 
occur,  but  quantitatively,  scattering  by  structure  with  dimensions  on 
the  order  of  optical  wavelengths  or  somewhat  larger  must  dominate  over 
Rayleigh  scattering."  (2:14). 

For  wavelengths  around  Ifjjn,  Rayleigh  scattering  is  not  likely  to 
occur  in  fatty  tissues  which  are  comparitively  homogeneous  on  a  scale  of 
micrometers.  This  scattering  would  be  due  to  fat  molecules  and  be 
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proportional  to  A  and  to  V  where  V  is  the  volume  of  a  molecule. 

This  scattering  should,  therefore,  be  weak.  Further,  since  breast 
tissue  in  younger  persons  tends  to  be  more  fibrous  than  fatty  in  nature, 
the  contribution  to  scattering  on  the  part  of  fatty  tissues  can  be 
considered  small  in  transillumination  applications. 

The  situation  is  quite  different,  however,  when  the  scattering 
properties  of  collagen  fibers  are  considered.  With  regard  to  the 
transmission  properties  of  these  fibers  in  the  dermis,  Parrish  writes: 

"The  data  indicates  that  the  Beer-Lambert  relation  is  invalid  for 
dermis,  and  that  transmittance  is  both  higher  and  more  forward 
directed  for  longer  wavelengths  over  the  region  between  0.5  and 
1.23^m.  These  observations  suggest  that  scattering  is  of  major 
importance  in  the  dermis."  (2:17) 

In  discussing  Findlay's  study  of  pig  dermis,  it  was  pointed  out  that 
collagenous  material  exhibited  little  absorption  over  the  entire  visible 
spectrum.  Thus,  from  the  above  remarks,  it  may  be  concluded  first  of 
all,  that  collagen  fibers  represent  a  primarily  scattering  medium  rather 


than  an  absorptive  one.  It  should  be  noted  that,  on  a  particle  scale, 
collagen  fibers  exhibit  an  identity  of  structure  not  only  throughout  the 
body  of  an  individual  but  across  species  as  well  (see  Figures  7  and  8). 
Thus,  it  is  considered  here  that  scattering  by  collagen  in  the  breast 
is,  to  a  very  good  approximation,  comparable  to  that  found  in  dermal 
collagen. 

So  far,  a  discussion  of  the  optical  properties  of  fibrous  breast 
tissue  has  led  to  the  conclusion  that  this  material  is  highly  scattering 
and  contibutes  very  little  to  absorption.  It  remains  to  both  quantify 
this  conclusion  (in  so  far  as  it  is  possible  to  do  so  given  the 
comparative  lack  of  data  in  this  area)  and  to  characterize  the  angular 
distribution  of  radiation  scattered  by  collagen  fibers.  The  problem  of 
radiation  distribution  will  be  considered  first  to  be  followed  by  a 
discussion  of  the  approximate  dimensions  of  the  scatterers. 

It  is  possible,  based  on  the  data  presented  in  Figure  9,  to  place 
some  bounds  on  the  distribution  of  scattered  radiation  from  collagen 
fibers.  The  data  in  Figure  9  graph  A  for  the  thin  sample  of  dermis  of 
thickness  O.U3fJn\  shows  that  for  longer  wavelengths  the  scattering 
pattern  of  the  tissue  as  a  whole  is  concentrated  in  the  forward 
direction.  Compared  to  the  thick  sample  of  dermis,  the  thin  sample 
reflects  much  more  clearly  the  scattering  properties  of  individual 
microscopic  scattering  centers.  This  is  due  to  the  fact  that  the  light 
propagating  through  the  thick  sample  is  multiply  scattered  many  times 
over  effectively  masking  or  averaging  out  the  contributions  of 
individual  scatterers.  The  data  previously  presented  in  Figure  5  giving 
percent  remittance  (R)  and  transmittance  (T)  for  a  20CfJn  thick  section 
of  dermis  is  reproduced  in  Figure  9  as  grapah  B.  Graph  B  of  Figure  9 
gives  absolute  values  of  remittance  and  transmittance  while  graph  A 


Figure  7.  Rat  tail  tendon  collagen  fibril  (39:16). 
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Figure  9.  Graph  A:  Scattering  Angle  Dependence  on  Wavelength 
(17:259).  Upper  curves  correspond  to  thin  sample;  lower  curves  to 
thick  sample.  Graph  B:  Spectral  transmittance  and  remittance  of  a 
200rfn  thick  section  of  human  dermis.  Data  is  given  in  percent 
remittance  (R)  or  transmittance  (T). 


represents  percentage  of  total  output  scattered  in  a  direction  at  the 
output  plane.  Of  course,  the  data  in  graphs  A  and  B  of  Figure  9  must  be 
consistent  with  one  another  and  an  adequate  theory  of  scattering  in  the 
breast  must  account  for  both  sets  of  data. 


An  assumption  that  the  microscopic  scattering  patterns  of  individual 
scattering  centers  are  forward  directed  can  be  shown  to  be  consistent 
with  all  the  data  presented  in  Figure  9.  First,  it  is  important  to  note 
that  this  data  was  collected  from  a  gross  sample  of  20<^m  in  thickness. 
The  transmittance  and  remittance  values,  therefore,  are  the  result  of 
scattering  contributions  from  a  great  many  scattering  centers.  The  fact 
that  there  are  many  scattering  centers  (consistent  with  the  intense 
scattering  known  to  characterize  meshes  of  collagen  fibers)  implies  that 
virtually  all  the  photons  are  scattered  during  propagation  through  the 
tissue.  Further,  graph  B  of  Figure  9  shows  that  at  approximately  1.2^m 
where  transmission  is  greatest,  the  ratio  of  transmitted  to  remitted 
light  is  about  four  to  one.  Therefore,  the  most  likely  explantion  for 
the  narrow  output  scattering  pattern  of  the  thin  sample  for  longer 
wavelengths  is  that  the  individual  scattering  centers  in  the  tissue 
scatter  light  in  a  primarily  forward  direction.  In  accordance  with  the 
general  principles  regarding  size  of  particles,  wavelength,  and 
scattering  outlined  on  pages  17  and  18,  this  forward  directed  scattering 
implies  that  the  scattering  centers  are,  on  the  average,  larger  than  a 
wavelength  of  about  yjfr u  The  work  of  Hardy,  et  al.  (17)  supports  this 
conclusion: 

For  thicknesses  as  great  as  1mm  scattering  is  maximal  at  all 
wavelengths  although  for  the  thinner  specimens  there  is  a  wavelength 
dependence  of  scattering.  This  dependence,  although  varying 
inversely  with  wavelength,  is  far  different  from  the  Rayleigh 
scattering.  This  is  to  be  expected  as  the  'particle  size'  of  the 
scattering  elements  of  the  skin  specimen  is  probably  not  less  than 
that  of  a  wavelength  of  light  in  the  spectral  range  studied  in 
these  experiments  (17:258). 
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At  first  glance,  however,  this  conclusion  regarding  the  scattered 
radiation  patterns  seems  to  be  at  variance  with  the  fact  that  collagen 
fibers,  though  they  may  be  a  few  micrometers  in  length,  have  an  average 
diameter  of  approximately  lOOnm.  This  fact  makes  each  fiber,  in  one 
dimension,  at  least  10  times  smaller  than  a  IfJ n  wavelength.  Since  the 
fibers  are  much  smaller  than  a  wavelength,  one  would  expect  the 
scattered  radiation  to  be  distributed  almost  equally  in  space  at  least 
about  the  axis  of  the  fiber.  To  conclude  that  the  scattered  radiation 
is  isotropic  in  the  axial  direction,  would,  however,  be  incorrect  as 
the  following  will  show. 

Collagen  fibers  derive  their  collective  strength  as  supportive 
structures  from  their  high  number  density  in  a  unit  volume.  Some  idea 
of  their  density  and  arrangement  may  be  gleaned  from  Figure  8  which 
shows  a  number  of  these  fibers  taken  from  a  human  specimen.  As 
particles  are  packed  closer  together,  the  scattering  properties  they 
exhibit  individually,  or  in  small  numbers,  change  dramatically 
(33:913-914).  For  dense  particles,  the  scattering,  in  fact,  decreases 
since,  the  denser  the  particles  get,  the  more  they  begin  to  resemble  a 
homogeneous  medium.  Under  these  circumstances,  the  product  of  the 
number  density  particles  of  O  and  the  scattering  cross  section  (7 
becomes p^CT  =  w(l  -  w)C^/Vg.  Here,  w  is  the  fraction  of  the 
total  volume  of  the  gross  medium  which  all  the  scattering  particles  take 
up  and  Vg  is  the  volume  of  a  single  particle.  The  factor  w(l  -  w)  is 
parabolic  with  a  maximum  value  at  w  *  0.5.  This  means  thatpnC^  has 
the  same  value  at  w  *  0.95  that  it  has  a  0.05.  However,  the  identity  of 
the  scatterers  themselves  may  change  from  that  of  the  particles  at 
w  *  0.05  to  the  interstitial  gaps  between  particles  at  w  1  0.95. 

Considering  Figure  8  and  drawing  upon  the  arguments  above,  one  is 


led  to  conclude  that  in  a  dense  mesh  of  collagen  fibers,  the  actual 
scattering  centers  are  probably  more  accurately  identified  with  the  gaps 
between  fibers  or  bundles  of  fibers  with  dimensions  comparable  to  or 
perhaps  larger  than  an  optical  wavelength.  This  estimate  in  size  of  the 
scattering  centers  should,  in  all  liklihood,  be  revised  upward  to  obtain 
a  large  particle  approximation  to  scattering  consistent  with  the 
transmittance  and  remittance  data  in  Figure  5.  Considering,  again,  the 
data  in  Figure  9,  it  may  be  felt  that  the  angular  spread  of  radiation 
for  shorter  wavelengths  (such  as  0. 55(^0  is  inconsistent  with  large 
particle  scattering  since  these  short  wavelengths  are  smaller  than  the 
dimensions  of  gaps  between  fibers.  A  possible  explanation  for  the 
phenomenon  may  be  that  the  shorter  wavelengths  are  exciting  a  different 
scattering  species  (perhaps  the  fibers  themselves). 

To  close  this  section,  a  final  comment  regarding  optimal 
transmittance  through  dermal  tissue  is  called  for.  The  reader  will 
recall  from  Figure  5  that  transmittance  is  greatest  at  a  wavelength  of 
about  1.2^Cin.  This  suggests  that,  if  internal  fibrous  tissue  of  the 
breast  is  similar  to  dermal  collagen  in  structure,  then  a  relatively 
good  value  of  transmittance  may  be  obtained  in  breast  transillumination 
at  1.2^/n.  At  this  wavelength  absorption  in  the  fibrous  tissue  is  low 
(about  20  percent)  and,  therefore,  light  contrast  due  to  differential 
absorption  between  cancers  and  surrounding  tissues  should  be 
comparitively  good.  The  suggestion  to  detect  light  of  1.2^Jn 
is  made  in  view  of  the  fact  that  detection  in  transillumination  is  now 
performed  using  silicon  detectors  whose  responses  extend  to  about  1.^2m 
but  no  further.  Therefore,  the  transmission  data  so  far  collected  in 
transillunination  does  not  include  transmittances  for  wavelengths 
beyond  l.^m. 


Although  transmittance  maxima  also  occur  at  1.7  and  2.2fjn, 
wavelengths  in  these  regions  suffer  greater  absorption  and  it  can  be 
expected  that  contrast  in  transillumination  at  these  wavelengths  would 


decrease  somewhat  from  that  at 


1.2^n. 


Light  Scattering  in  Whole  Blood. 

Whole  blood  consists  of  red  blood  cells  and  blood  plasma.  In  normal 
blood  tissue,  the  red  blood  cells  account  for  approximately  40  percent 
of  the  total  blood  volume.  Thus,  60  percent  of  blood  is  plasma  which  is 
almost  entirely  made  up  of  water. 

Considered  in  isolation,  red  blood  cells  which  are  about  ^jm  in 
diameter  and  ^jn  in  thickness,  are  highly  scattering  in  a  forward 
direction  at  near  infrared  wavelengths  (19:63,66).  However,  because  the 
cells  are  surrounded  by  plasma,  radiation  spanning  the  red  and  near 
infrared  is  absorbed  in  the  plasma  either  before  it  reaches  the  cells  or 
after  scattering  from  them.  This  accounts  for  the  shadows  cast  in 
trans illumination  by  superficial  blood  vessels  and  by  the  mass  of 
proliferated  blood  vessels  which  surround  the  highest  percentage  of 
tumors  in  the  breast.  Consequently,  whole  blood  in  breast  tissue  is 
considered  primarily  an  absorptive  medium  vice  a  scattering  one  in 
trans i 1 lumina t ion . 

Taking  together,  then,  the  data  and  arguments  presented  above,  one 
may  obtain  a  reasonable  model  of  light  interaction  with  breast  tissues 
at  red  and  near  infrared  wavelengths.  The  features  of  this  model  would 
include 

(1)  Low  absorption  in  tissues  surrounding  blood  vessels  and 
cancers , 

(2)  High  absorption  in  blood  and  blood  vessels  massed  about 


cancers , 


(3)  Low  scattering  in  fatty  and  whole  blood  tissues,  and 

(4)  Intense  highly  forward  directed  scattering  in  the  intralobular 
connective  tissue  made  up  of  collagen  fibers. 

This  model  will  be  assuned  in  Chapter  II  in  the  development  of  the 
mathematical  theory  of  multiple  scattering  as  it  applies  to 
transillunination.  There  the  reader  will  find  the  propagation  of  near 
infrared  wavelengths  through  tissues  surrounding  breast  tumors  to  be 
described  by  mathematical  functions  which  provide  for  forward  directed 
scattering,  relatively  large  scattering  cross  section  0^,  and  low 
absorption. 


II.  Scattering  Theory 

As  will  become  clear  in  the  discussion  to  follow,  a  central  goal  of 
the  mathematical  development  of  scattering  theory  as  it  applies  to 
breast  transillumination  must  be  to  calculate  the  correlation  function 
(or  mutual  coherence  function)  for  a  scattering  medium.  Ihe  correlation 
function  is  an  important  quantity  in  the  calculation  of  both  the  point 
spread  function  and  the  limitation  of  resolution  of  an  imaging  system 
in  the  presence  of  a  random  medium.  It  will  become  clear  that  the 
correlation  function  is  the  key  to  a  quantitative  determination  of  the 
quality  of  imaging  that  can  be  expected.  From  a  knowledge  of  the  point 
spread  function,  one  can  estimate  the  efficacy  of  the  techniques  of 
image  processing  such  as  phase  retrieval  which  employ  extrapolation 
from  the  Fourier  modulus. 

The  mathematical  development  to  follow  proceeds  along  two  differing 

lines  of  thought  corresponding  to  the  analytic  and  transport  theories  of 

scattering  in  a  random  medium.  The  analytic  theory  (also  called 

multiple  scattering  theory)  begins  with  Maxwell's  equations  or  the  wave 

equation  and  develops  differential  or  integral  equations  for  statistical 

quantities  of  interest.  Ihe  transport  theory,  on  the  other  hand,  is 

not  based  on  an  electromagnetic  wave  approach  but  deals  directly  with 

the  movement  of  energy  through  a  medium  containing  scattering  and 

absorbing  particles.  Since  the  transport  theory  does  not  develop  from 

electromagnetic  wave  propagation,  it  does  not  directly  include 

diffraction  and  interference  effects — these  effects  enter  indirectly 

through  the  scattering  cross  section  (Cf )  and  the  absorption  cross 

section  (C7)  (CT  and^T  being  measurable  quantities)  (19:147). 
a  s  a 

Although  the  analytic  and  transport  theories  of  scattering  take 


quite  different  points  of  view,  it  can  be  shown  that  the  two  theories 
are  connected.  Specifically,  the  correlation  function  of  analytic 
theory  and  the  specific  intensity  of  transport  theory  are  a  Fourier 
transform  pair  (20:275).  Thus,  although  the  transport  theory  does 
not  contain  information  as  to  fields,  it  is  possible  to  calculate 
the  correlation  of  fields  via  a  Fourier  transform.  This  concept  will 
be  exploited  when  the  correlation  function  for  a  slab  of  tissue  will 
be  derived  based  on  a  knowledge  of  the  specific  intensity. 

Twersky ' s  multiple  scattering  theory  is  developed  in  Appendix  A 
where  the  integral  equation  for  the  correlation  function  is 
obtained.  This  theory  was  chosen  as  it  has  been  successfully 
applied  by  investigators  to  scattering  problems  in  biological 
media  (30).  Twersky's  correlation  function  is  defined  and  its  relation 
to  the  specific  intensity  of  transport  theory  through  a  Fourier 
transform  is  shown. 

A  solution  of  the  transport  equation  as  tailored  to  the  case  of  the 
plane-parallel  slab  of  tissue  is  obtained.  The  slab  configuration  was 
chosen  as  it  corresponds  well  to  the  clinical  set-up  for 
transillumination  of  the  breast  where  the  breast  is  itself  compressed 
into  the  form  of  a  plane-parallel  mediun  (precisely  as  is  done  in  X-ray 
mammography). 

This  chapter  presents  four  topics:  (1)  the  correlation  function  for 
analytic  theory,  (2)  the  connection  between  the  analytic  and  transport 
theories,  (3)  the  limitation  on  image  resolution  due  to  a  scattering 
mediun,  and  (4)  the  signal  to  noise  ratio  relevant  to  transillumination. 
It  will  be  shown  that  according  to  scattering  theory  embodied  in  the 
four  topics  above,  the  image  obtained  via  transillumination  cannot  be 
improved  without  first  reducing  the  scattering  by  some  physical  means 


before  image  processing  is  applied.  Ibis  conclusion  will  lead  in 
Chapter  III  to  a  consideration  of  time  gating  and  spatial  filtering-  a 
method  of  imaging  which  attempts  to  reduce  the  effective  scattering 
(and  thus  reduce  noise)  by  temporal  and  spatial  discrimination  of  photon 
paths. 

The  Correlation  Function  From  Analytic  Theory 

The  correlation  function  I^r^,  r^)  of  analytic  theory  for  the 
scattered  field  at  points  and  is  presented  here  as  it  will  be  used 
later  to  calculate  the  limit  of  resolution  one  can  expect  to  obtain  when 
imaging  through  a  random  medium.  Ibis  limitation  of.  resolution  is 
important  to  the  trans illumination  problem  as  it  will  define  the 
theoretical  limits  of  the  efficacy  of  transillumination  and  will  provide 
quantitative  methods  for  determining  these  limits. 

Twersky's  integral  equation  for  the  correlation  of  the  field \pa  at  a 
point  ra  and ^  at  a  point  in  a  random  medium  is  (35:99) 


<^v> = f ></■> 
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vg  and  vg  are  called  multiple  scattering  operators.  Tbeir  function 

is  to  convert  a  wave  incident  on  a  scatterer  at  a  location  r  into  a 

s 

scattered  wave  which  is  then  subsequently  observed  at  location  r  (in 

3 

3  b 

the  case  of  vg  )  or  (in  the  case  of  vg  ).  Tbe  function  P(rg) 
represents  the  number  density  of  scatterers  per  unit  volume  as  a 
function  of  position.  In  equation  (2)  the  symbol  uga  represents  the 


scattering  of  a  wave  incident  on  a  scatterer  at  location  rg  and 

propagating  through  free  space  to  a  point  r  .  Thus,  u  is  called  the 

3  b 

free  space  scattering  operator.  The  form  of  the  operators,  v  ,  v  , 

s  s 

and  ug  are  derived  in  detail  in  Appendix  B.  The  brackets  <  >  in 
equation  (1)  indicate  the  usual  time  average  value  and  equation  (1)  is 
integrated  over  all  scatterers  in  the  medium.  In  order  to  complete  the 
calculation  of  the  correlation  function  given  by  equation  (1),  it  would 
be  necessary  to  find  expressions  for  the  quantity  >  and 

<|^s|2>  using  the  analytic  theory.  But  complete  solutions  for  the 
Twersky  integral  equation  for  the  correlation  function  have  not  appeared 
in  the  literature.  However,  as  adumbrated  at  the  beginning  of  this 
chapter,  a  methodology  does  exist  for  the  complete  calculation  of  the 
correlation  function  based  on  the  fact  that  the  correlation  function  of 
analytic  theory  and  the  specific  intensity  of  the  transport  theory  are  a 
Fourier  transform  pair.  It  remains  now  to  justify  the  assertion  that 
this  Fourier  transform  connection  exists  between  the  analytic  and 
transport  theories.  Once  this  connection  is  established,  it  will  be 
possible  to  depart  from  the  analytic  theory  and  calculate  the 
correlation  of  fields  from  the  specific  intensity.  The  reader  will 
find  a  description  of  the  transport  theory  and  specific  intensity  in 
Appendix  C. 

The  Connection  Between  the  Analytic  and  Transport  Theories 

Although  solutions  for  Twersky' s  correlation  function  based  on 
analytic  theory  have  not  yet  appeared,  various  solutions  for  the 
specific  intensity  in  radiation  transport  theory  do  exist.  The  reader 
is  referred  to  Appendix  C  for  a  description  and  discussion  of  transport 
theory.  In  this  section,  the  connection  between  the  correlation 
function  and  specific  intensity  is  described.  With  this  knowledge  in 


hand,  it  will  be  possible  to  obtain  information  regarding  correlation 
of  fields  despite  the  fact  that  radiation  transport  theory  is  not  based 
on  electromagnetic  field  equations.  In  later  developments,  it  will  be 
seen  that  the  correlation  function  enters  as  an  essential  quantity  in 
the  calculation  of  the  limit  of  resolution  of  an  imaging  system.  Thus, 
the  transport  theory  becomes  a  powerful  tool  in  the  analysis  of  imaging 

of  the  breast  by  trans illumination.  The  specific  intensity  I(r,s)  is 

-2  -1  -1 

measured  in  Wm  sr  Hz  and  is  the  average  power  flux  density  within  a 
unit  frequency  band  centered  at  frequency  and  radiating  in  the 
direction  s  from  a  point  r  on  a  radiating  surface.  See  Figure  10. 


Figure  10.  Flux  through  da  on  a  surface  A  (19:150). 

The  correlation  function  presented  in  equation  (1)  and  the  specific 


intensity  are  related  by  (20:275) 

=  ITtojOb)}  =r (?,$>  =  f  I(f,§)exp(iKr§.rd)(jn  (3) 

where  r  =  l/2(ra  +  r^)  and  r^  =  (r  -  r^),  and  Kr  is  the  real  part  of 
the  complex  wave  number  for  a  scattering  medium  (see  Appendix  B  for  a 
description  of  K). 

It  can  be  shown  that,  assuming  the  validity  of  equation  (3),  the 
integral  equation  for  the  average  intensity  in  transport  theory  can  be 
derived  (see  Appendix  D).  Since  this  is  the  case,  it  follows  that  the 


correlation  function  of  analytic  theory  and  the  specific  intensity  of 
transport  theory  are  a  Fourier  transform  pair  precisely  as  stated  in 
equation  (3). 

Limitation  on  Image  Resolution  Due  to  a  Scattering  Medium 

The  reader  will  recall  from  Chapter  I  that,  for  the  purposes  of 
transillumination,  the  breast  is  compressed  into  a  plane-parallel 
optical  medium.  In  this  section,  the  limit  of  resolution  of  an  imaging 
system  in  the  presence  of  a  plane-parallel  scattering  mediun  is 
considered.  The  reason  for  calculating  such  a  resolution  limit  is  that 
it  provides  a  measure  for  assessing  the  efficacy  of  imaging  via 
trans illumination.  How  this  is  the  case  will  become  clear  as  the 
discussion  proceeds.  Using  the  Fourier  transform  relation  between  the 
correlation  function  and  the  specific  intensity,  the  correlation 
function  for  a  plane  wave  (in  the  large  particle  scattering 
approximation)  incident  upon  a  scattering  medium  is  presented.  Figure 
11  depicts  the  physical  set  up  for  the  calculation  of  the  limit  of 
resolution. 


II  - 

« 

o  0<0  a 

iWy*) 

M 

(juy) 

X 

II 

«  0  * 

s 

•4 

« 

O  a  ° 

0  C» 

6  °  0 

1 ,  r 

*  £> 

1 

1 

1 

Figure  11.  A  monochromatic  plane  wavej//.  propagating  through  a 
random  distribution  of  scatterers  (20:301).  1The  resultant  scattered 
wave  is  focused  cuto  the  x,y  plane  (Fourier  transform  plane).  The  lens 
has  radius  a  and  focal  length  f. 


The  assumption,  based  on  the  optical  model  presented  in  Chapter  I, 


that  scattering  is  due  primarily  to  particles  larger  than  a  wavelength 
in  the  near  infrared  implies  that  the  scattering  pattern  of  an 
individual  particle  is  forward  directed.  A  mathematical  formulation 
of  a  forward  scattering  pattern  may  be  obtained  by  means  of  a  phase 
function  p(s)  which  describes  the  amount  of  power  which  is  scattered 
into  a  direction  s  by  a  scatterer.  A  forward-directed  scattering 
pattern  may  be  approximated  by  giving  p(s)  a  Gaussian  form: 


pcsj  -  4apW0 exp (-OpS1) 


where,  for  large  scatterers  Ot^  -  Z.UL  where  D  is  the 

particle  diameter  and  X  is  the  wavelength  of  incident  radiation.  W q  is 
called  the  albedo  of  the  scatterer  and  is  defined  as 


V  csAas  +  cra)  =  «{/<t 


where  CT  is  the  scattering  cross  section  of  a  particle,  O'  is  the 
absorption  cross  section,  and  Oj.  (total  cross  section)  is  equal  to 

0+0. 

s  a 

In  appendices  E  and  F,  it  is  shown  that  the  correlation  function 
^(z,q)  for  a  plane-parallel  medium  with  plane  wave  incidence  is 

r(H,q1  =  I„exp  |-/0£^e[  1  -  W0exp(-  . 

m 


where  Pn  is  the  number  density  of  scatterers  per  unit  volume.  The 
quantity  q  =  K^yQj '  where  Kr  is,  again,  the  real  part  of  the  complex  wave 
nunber  for  a  scattering  median  and/^'  =  ^  '  -  @2 '  where  1  and  ' 
are  points  in  the  x',  y'  plane  as  picture  in  Figure  11.  Thus,  I~(z,q) 
represents  the  correlation  of  fields  originating  at  points  and 


as  measured  a  distance  z  in  the  horizontal  direction  in  Figure  11. 


The  intensity  in  the  focal  plane  of  the  lens  pictured  in  Figure  11 
is  equal  to  the  point  spread  function  P ^(f))  for  the  system  depicted 
since,  in  this  case,  the  image  is  that  of  a  point  source  at  infinity. 

In  physical  terms,  the  point  spread  function  describes  the  effect  the 
scatterers  have  on  the  resultant  image.  Heuristically,  the  point 
spread  function  states  that  the  point  source  is  "spread  out"  in  the 
focal  plane.  The  reader  will  note  the  conspicuous  absence  of  a  cancer 
in  the  medium  in  Figure  11.  The  reason  for  this  is  that  the  goal  of 
calculating  the  resolution  limit  is  to  gain  an  understanding  of  the 
effect  of  the  scatterers  on  a  propagating  wave.  The  image  that  results 
from  imaging  an  opaque  object  (in  this  case  a  cancer  surrounded  by  blood 
vessels)  through  a  non-scattering  mediun  is  well  understood.  Thus, 
it  is  only  necessary  to  image  the  scatterers  alone  to  obtain  the 
required  understanding  of  the  effect  of  scattering  on  cancer  detection 
by  transillunination. 

The  point  spread  function  P^(0)  or  intensity  in  the  focal  plane  of 
Figure  11  can  be  shown  to  be  (see  Appendix  G) 


X(W)nz,3')K(p/) 


(7) 


where  k  =  f  =  focal  length  of  a  lens,  a  *  radius  of  lens,  Kp’^) 

is  the  pupil  function,  and  ~  quantity p  is  a  point  x,y 

2  2  1/2 

at  a  distance  (x  +  y )  '  from  the  horizontal  axis  (see  Figure  11). 

In  Appendix  G  it  is  shown  that  for  an  optically  thick  medium,  the 
point  spread  function  P^p)  may  be  expressed  as  the  sum  of  a  coherent 
(average)  part  and  an  incoherent  (fluctuating)  part: 


(8) 


Pf  (P)  =  PC(P)  +  \(p) 


where  P c<p)  is  the  coherent  portion  of  P^(p)  and  P ^{f))  is  the  incoherent 
portion.  Also  shown  in  Appendix  G  are  expression  for  P ^(p)  and  P 


?C(P)~  l0exp(-T)  (a/p)2  Jj2 (*pa/f ) 


?i(P)  =  I^a/Pifexp  [-r(i-  Vl0)  -  (P/P^  ] 


whereT; optical  thickness  Again,  Pn  is  the  nunber 

density  of  scatterers  per  unit  volume,  CT  and  CT  are  the  scattering  and 

S  3 

absorption  cross  sections,  and 


P*  =  (Kr£/x)2(TWc/ap) 


The  quantity  z  in  the  expression  for  Tis,  for  the  transillumination 


problem,  the  physical  thickness  of  a  breast  compressed  for 


transillumination.  The  relationship  between  Pc(fi)  and  P ^(p)  is 
shown  in  Figure  12  where  P^(p)  is  superimposed  over 


Figure  12.  P  (0)  is  the  Airy  pattern  created  by  the  imaging  of  the 
coherently  propagating  portion  of  the  intensity  (20:305).  The  incoherent 
intensity  P .(p)  is  spread  out  over  the  focal  plane  and  increases  with 
increasing  optical  thickness  T. 

The  coherent  intensity  is  the  image  bearing  component  of  the  total 
intensity.  If  P^(0)  is  less  than  P c(p)  around  the  pointp=  0,  it  is 
possible  to  obtain  an  image  of  an  object  surrounded  by  a  scattering 
medium.  However,  if  P ^(ft)  is  greater  than  P C<P),  the  image  is  "washed 
out"  by  the  incoherent  intensity  (20:305).  Note,  however,  that  even  in 
cases  where  P ^(p)  <  Pc<p),  if  P^p)  is  significant,  the  contrast  in 
the  image  may  be  quite  poor.  Under  these  circumstances  some  resolution 
may  be  possible  but  the  image  will  still  be  difficult  to  see. 

The  Signal  to  Noise  Ratio 

Based  on  the  development  in  the  previous  section  it  is  possible  to 
define  a  signal  to  noise  ratio  in  terms  of  the  coherent  intensity  P c<p) 
and  the  incoherent  intensity  P $))-  It  should  be  understood  at  the 
outset  that  the  noise  (which  will  be  defined  as  the  incoherent 
intensity)  is  not  additive.  In  fact,  the  incoherent  noise  is  created 
directly  from  the  coherently  propagating  intensity  as  it  is  scattered  by 
the  medium.  Thus,  as  the  incoherent  intensity  grows,  the  coherent 
intensity  diminishes. 


The  signal  to  noise  ratio  will  be  defined  as  the  ratio  of  the 
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coherent  to  incoherent  intensities: 


PC(P)  ..  ..  s 

PjCO  6  n 


(12) 


To  simplify  the  calculation,  consider  the  signal  to  noise  ratio  at  the 


point  P=  0: 

s 
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Expanding  the  numerator  using  the  formula 


J^x)  a  x/z  -  X/424  +  x%z4z6  X  4  £  8 
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and  then  recalling  that 


,.2  a  (Kr  fYk2-)  (rWo/a  ) 


(U) 
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the  signal  to  noise  ratio  becomes 


(Kra/zf  (rW0/ap)exp  (-rW0) 


(15) 


It  can  readily  be  seen  that  the  signal  to  noise  ratio  is  strongly  a 
function  of  T  —  P^ as  S/N  falls  exponentially  with 
increasing  optical  thickness. 

As  an  example  of  a  calculation  of  S/N  in  a  biological  medium,  the 
signal  to  noise  ratio  for  a  dilute  blood  medium  will  be  obtained.  The 
physical  parameters  employed  are  the  same  as  those  used  by 
Maarek  et  al.  (24)  in  a  Monte  Carlo  simulation  of  photon  propagation  in 
a  tissue  medium  (see  the  section  entitled  '*Time  Gating  the  Output"  in 
Chapter  V).  This  allows  comparison  of  calculation  techniques. 

Before  the  calculation  is  attempted,  however,  a  brief  digression  to 
discuss  the  relationship  of  the  scattering  cross  section  to  particle 
density  is  required.  As  the  number  density  of  scatterers  increases,  the 
particles  are  packed  closely  together.  As  the  ratio  of  total  volune  of 
scatterers  to  total  volume  of  medium  approaches  1,  the  particles  are 
packed  together  in  such  a  manner  as  to  create  a  homogeneous  medium  of 
particles  in  which  no  scattering  occurs.  Thus,  based  on  the  argument 
above,  one  is  led  to  the  following  approximation  for  the  factor /^n^: 

Pncr  =[v»(i-W)<%]/ve  as) 

where  w  is  the  ratio  of  the  total  volume  of  particles  to  the  total 
volume  of  the  medium  and  Vg  is  the  volume  of  a  single  particle. 

Equation  (16)  was  obtained  by  Twersky  via  a  more  mathematically 
rigorous  approach  (33:913).  Incidentally,  these  considerations  throw 


some  light  on  why  transillumination  can,  in  some  instances,  give  better 
results  for  dense  tissues  than  can  x-ray  mammography. 

As  in  the  study  done  by  Maarek  et  al.  (24)  the  following  calculation 
assumes  a  hematocrit  H  (the  ratio  of  blood  cell  content  to  whole  blood 
where  vdiole  blood  consists  of  blood  cells  and  plasma)  of  0.05.  Thus,  in 
this  case,  the  quantity  H  =  w.  In  addition  to  H,  the  equations  listed 
below  were  substituted  into  equation  (15)  to  calculate  the  signal  to 
noise  ratio: 

W0  =  <J/St  =■  + 

Kr  =  K  +  H/(  1+  H}  k[  n‘/nc  -  1] , 

ap  --  2.QG  (D/Xf , 

=  (e)+/g%z 

and 

Vc  =rrrzh 

where  r  is  the  radius  of  a  blood  cell  and  h  is  its  thickness.  The 
equation  for  the  real  part  of  K  (Kf)  is  taken  from  Twersky  (32:168) 
and  the  ratio  n^/nQ  is  the  ratio  of  the  index  of  refraction  of  a  blood 
cell  to  that  of  water.  As  will  be  seen,  the  accuracy  of  the  second 
term  used  in  Kr  are  not  critical  as  the  signal  to  noise  ratio  is  so 
strongly  dependent  on  T  . 


The  numerical  values  chosen  were 
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z  =  tissue  thickness  =  4cm 

The  values  for  G  and  G  are  those  supplied  by  Ishimaru  (19:66)  while 
blood  cell  dimensions  (D,h)  are  taken  from  Johnson  and  Guy  (21).  The 
wavelength  (\)  is  consistent  with  Ishimaru's  data  on  scattering  and 
absorption  cross  sections. 

Substitution  of  these  numerical  values  into  equation  (11)  gives  the 
following  result  for  the  signal  to  noise  ratio: 


(Z.21  xl01°)cxp(-74G.lG) 

2.21x10‘334-O 


(17) 


Equation  (17)  indicates  that  the  coherent  intensity  has  been  completely 
turned  into  incoherent  intensity  by  the  scattering  process.  It  might 
now  be  instructive  to  inquire  at  what  point  (in  terms  of  optical 
thickness)  one  might  expect  to  obtain  a  coherent  image.  This  question 

pc  p) 

may  be  answered  by  setting  the  ratio  -  equal  to  unity  and  then 

Pi<£)  =  0 

solving  for  the  optical  thickness  T. 

Using  equation  (15)  one  obtains 
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vc(p) 

Pi  <p) 


I  p:0 


KtfafrWp 


expC-rW,) 


=  1 


which  implies  that 


(18) 


The  quantity  Ct-  is  proportional  to  (D/\)  and  K^ssk.  Therefore,  the 

*  2 

left  side  of  equation  (18)  is  proportional  to  (a/D)  —  a  quantity  which 

is  likely  to  be  quite  large  in  the  case  of  trans illumination.  It  may  be 

concluded,  then,  that  equation  (18)  is  true  only  for  large 

(i.e.,  either  strong  scattering  and/or  large  thickness  of  mediun). 

Letting  (Kra)X_  V)  and  tIaL  -  X  ,  the  following  mathematical  steps 
40p  “  V 


establish  an  equation  for  optical  thickness  wherein  one  may  expect  to 
obtain  a  reasonably  good  image  in  a  scattering  medium.  First,  following 
the  definitions  of  B  and  x  given  above,  equation  (18)  implies  that 

B  =  expcx)/x 


Therefore, 

inb  +  In.  (in  S')  ~  InB  +  Itl  (x- lnx) 

-  ln[B(*-l^x)] 

-  Tn[Bx-Blnx] 

Z  lnCcxpCx^-^KpCxVx)!^^] 

r  WxTl 

If  x  >>  lnx,  then 


1tiB+  lne*  -  X  (19) 

where  and  x=rW0 . 
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Now  if  the  diameter  D  of  a  scatterer  is  large  compared  to  A,  then&p  = 

2.66  (D/^)^  (19:122).  Substituting  =  2.66  (D/A)^  into  the  equation  for  B 
yields 


■B  =  (Kr^T  =  3.71  Ca/D)2 
4ap 

Letting  a  =  2cm  and  D  =  j^m  (as  in  the  example  calculation  for  S/N)  one 
obtains 

"B  =  3.05  x  107 


Substituting  this  value  of  B  into  equation  (19)  one  gets  that 

In  (3.03  x107)+1h(1ti  3.03x107)=  rW0  ( 

For  a  medium  that  is  not  highly  absorbing  (as  is  the  case  with  breast 
tissues  other  than  blood)  WqSI  andTW^sT.  Therefore,  from  equation 
(21), T-  20.08  and 


=2aos 


(22) 


Then  using  the  fact  that  yOCF=  ^ ~  anc^  the  same  Parameter  values 
used  before  one  obtains  for  z  (thickness) 

Z  -  1.019mm 


This  value  for  z,  then,  represents  the  upper  limit  on  the  thickness  of 
tissue  for  which  one  would  expect  to  be  able  to  obtain  a  reasonable 
coherent  image,  i.e.,  an  image  which  would  still  contain  details  of  the 
object  profile. 


Interpretation  of  the  Signal  to  Noise  Ratio. 

The  calculation  of  S/N  shows  that,  for  the  parameters  chosen,  the 
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coherent  signal  Pc(p)  is  (for  all  intents  and  purposes)  non-existent. 

_334 

The  exponential  factor  (10  )  may  seem  surprising  since  in  many 

cases  of  noisy  systems,  although  the  signal  to  noise  ratio  may  be  quite 
low,  equation  (17)  indicates  an  extremely  weak  signal,  indeed  (if  it 
may  be  said  to  be  there  at  all).  However,  if  one  considers  the  fact 
that  the  noise  is  not  additive  but  derives  from  the  scattering  of  the 
coherent  signal  itself,  it  is  easier  to  see  how  the  coherent  intensity 
should  be  so  overcome  by  the  incoherent  intensity  that  it  "feeds".  If 
the  scattering  cross  section  is  at  all  significant,  it  does  not  take 
long  as  the  wave  propagation  through  the  medium  for  the  coherent 
portion  of  the  intensity  to  be  removed  almost  completely. 

Given  the  analysis  presented  in  this  chapter,  it  does  not  appear 
that  the  optical  processing  techniques  aimed  at  retrieving  or  enhancing 
the  coherent  intensity  in  order  to  improve  imaging  and  are  likely  to 
succeed.  Such  methods  as  matched  filtering  or  simple  averaging  of 
outputs  to  eliminate  or  reduce  noise  assume  that  the  signal  is  present 
in  the  noise  and  is,  therefore,  retrievable,  at  least  in  principle.  In 
the  case  of  breast  transillumination,  however,  the  coherent  or  average 
signal  is  no  longer  present  in  the  total  wave.  Therefore,  a  matched 
filter  has  virtually  nothing  to  match.  On  the  other  hand,  although 
averaging  outputs  (at,  say,  different  angles)  could  theoretically  reduce 
noise,  the  result  would  not  be  an  enhanced  coherent  signal  (since,  in 
actuality,  there  virtually  is  none).  Of  course,  the  parameters  chosen 
here  for  the  calculation  of  the  signal  to  noise  ratio  do  not  necessarily 
correspond  exactly  to  those  of  breast  tissue.  The  values  of  the 
parameters  were  chosen  largely  as  they  were  simply  because  the 
scattering  and  absorption  cross  sections,  density  and  other  physical 
parameters  for  the  breast  were  unavailable.  Nevertheless,  it  is 
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reasonable  to  assume  that  breast  tissue  would  not  exhibit  properties  so 
different  from  those  considered  in  this  chapter  that  the  conclusion 
implied  by  equation  (17)  would  be  changed  in  substance.  The  reader  will 
also  find  that  the  calculation  of  the  total  intensity  to  be  accomplished 
in  Chapter  III  will  show  that  the  parameters  used  here  yield  results 
fairly  close  to  those  obtained  under  actual  transillumination. 

Further,  one  should  consider  that  the  assumption  of  a  primarily 
forward  scattering  pattern  for  single  scatterers  constituted  a  "best 
possible"  scattering  medium  for  the  retrieval  of  the  coherent  intensity. 
If,  by  chance,  the  single  scatterers  in  the  breast  had  scattered 
radiation  patterns  broader  that  that  indicated  by  a  Gaussian 
approximation  to  the  phase  function,  the  signal  to  noise  ratio 
calculated  under  that  circumstance  could  only  be  worse  than  that  arrived 
at  in  the  analysis  presented.  That  is,  the  broader  the  scattering 
pattern  of  each  scatterer,  the  more  quickly  would  the  coherent  intensity 
be  spread  out  through  the  medium  and  the  more  quickly  would  the  coherent 
image  bearing  intensity  disappear. 

The  development  of  scattering  theory  presented  in  this  chapter  has 
led  to  the  conclusion  that  the  coherent  intensity  is  effectively 
destroyed  by  the  intense  scattering  process  in  a  biological  medium  of 
the  kind  considered  here. 

At  this  juncture,  it  is  important  to  formulate  the  problem  of 
transillumination  imaging  in  the  light  of  the  theoretical  results 
obtained.  In  particular,  it  is,  of  course,  experimentally  established 
that  images  of  tumors  are  observed  under  transillumination  despite  the 
fact  that  the  scattering  theory  presented  here  demonstrates  an  effective 
zero  value  of  coherent  image  bearing  intensity.  A  conclusion  that  may 
be  drawn  from  this  observation  is  that  the  imaging  occuring  in 


transillumination  is  not  obtained  by  means  of  the  coherently  propagating 
intensity.  Instead,  the  shadow  image  cast  by  a  tunor  must  be  the  result 
of  differences  in  incoherent  intensities.  Therefore,  to  improve 
transillumination  imaging  in  the  near  infrared,  one  should  first 
concentrate  on  improving  contrast  in  intensity  over  the  image  plane 
rather  than  attempting  to  resolve  the  details  of  the  object  of  interest 
(a  tumor)  imbedded  in  the  tissue. 

In  Appendix  H,  the  reader  will  find,  as  an  example  of  an  attempt 
to  recover  the  coherent  signal,  a  description  of  the  method  of  phase 
retrieval  as  it  might  be  applied  to  breast  transillumination.  The 
results  obtained  in  Appendix  H  show  that  the  phase  of  the  signal 
cannot  be  recovered  by  this  method.  Consequently,  phase  retrieval  does 
not  show  promise  of  improving  imaging  by  transillumination. 


III.  A  Method  for  Improving  Imaging  in  Transillunination 


In  this  chapter  a  method  for  improving  the  contrast  in 
transillunination  will  be  presented.  This  technique,  called  Time 
Gating  and  Spatial  Filtering  (TGSF),  concentrates  on  improving  contrast 
by  physically  reducing  the  effects  of  scattering.  The  method  does  not 
attempt  to  retrieve  a  coherent  signal  but  accentuates,  instead,  the 
contrast  in  incoherent  intensity  across  the  output  plane. 


Time  Gating  and  Spatial  Filtering  of  the  Transmitted  Radiation 


Recent  studies  (24)  based  on  a  Monte  Carlo  model  of  propagation  of 
light  through  thick  (20  mm  -  80  mm)  heterogeneous  biological  samples 
have  provided  encouraging  results  suggesting  that  spatial  and  time 
discrimination  of  output  radiation  can,  in  some  instances,  significantly 
improve  the  contrast  in  transillunination.  The  work  accomplished  thus 
far  has  not  been  directed  specifically  toward  the  problem  of  cancer 
detection  in  the  breast  but  the  method  used  is  certainly  applicable  to 
the  problem. 

Simulation  of  Laser  Tomoscopy 

This  section  is  devoted  to  a  description  of  the  Monte  Carlo 
simulation  of  photon  propagation  through  a  tissular  medium  performed  by 
Maarek,  et  al,  in  1986  (24).  In  this  study,  the  researchers  simulated 
the  illumination  of  a  tissular  slab  with  a  picosecond  laser  source. 
Inside  the  tissue  there  was  considered  to  be  a  blood  vessel  which  could 


be  moved  about  to  simulate  the  process  of  scanning  the  tissue  with  the 
source.  At  the  output  end  of  the  slab  was  a  simulated,  collimated, 
ultrafast  optical  Kerr  Shutter.  Thus,  at  the  shutter,  both  spatial  and 
temporal  discrimination  of  photons  could  be  performed  based  on 


scattering  angle  at  the  shutter  and  flight  times  of  photons. 


The  advantage  of  using  such  an  apparatus  can  be  understood  in  terms 
of  the  source  of  noise  in  the  imaging.  As  theory  indicates  (Chapter  II) 
and  experiment  demonstrates,  light  scattering  produces  a  reduction  in 
contrast  and  resolution  of  an  image.  In  terms  of  photon  paths,  the 
noise  at  the  detector  comes  from  photons  which  are  first  scattered  far 
from  the  axis  of  illumination  and  are  then  eventually  scattered  back 
into  the  field  of  view  of  the  detector.  These  photons  cause  light  areas 
to  appear  where,  imaging  in  the  absence  of  scattering,  there  might 
otherwise  be  dark  regions  (24:207).  Thus,  a  reduction  in  contrast  and 
resolution  occurs  in  the  image.  Now  those  photons  which  are  scattered 
far  from  the  axis  of  illumination  and  yet  reach  the  detector  travel 
distances  greater  (and  have  longer  flight  times)  than  those  photons 
which  tend  to  stay  nearer  the  illumination  axis.  Consequently,  a  time 
gate  at  the  output  slab  could  remove  noise  from  the  total  system  by 
"discarding"  those  photons  which  take  some  time  longer  than  a  pre-set 
acceptance  time  to  cross  the  tissular  slab.  The  time  gate  would  ensure 
that  only  those  photons  which  tend  to  travel  near  the  illumination  axis 
would  be  detected. 

Description  of  the  Simulated  Model 

The  tissular  model  which  Maarek,  et  al.  used  is  pictured  in  Figure 
13.  The  tissue  was  a  semi-infinite  slab  bounded  by  parallel  planes  at 
x  =  0  and  x  =  x  .  The  tissue  surrounding  the  blood  vessel  of  radius 
R  was  assimed  to  be  a  homogeneous  medium  with  optical  characteristics 
equivalent  to  that  of  dilute  blood  with  a  hematocrit  of  0.05  (indicating 
a  ratio  of  total  blood  cell  volune  to  blood  cells  and  plasma  of  0.05). 
The  blood  vessel  of  variable  radius  R  centered  at  coordinate  x  ,  y 
as  shown  in  Figure  13  contained  blood  with  a  hematocrit  of  0.50.  This 
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model  has  some  similarity  to  the  model  of  the  human  breast  described  in 
Chapter  I.  It  will  be  seen  later  in  the  description  of  the  results  of 
the  simulation  that  the  output  intensity  is,  in  terms  of  important 
qualitative  features,  very  similar  to  that  seen  in  transillumination  of 
the  breast.  In  both  the  surrounding  medium  and  the  blood  vessel  the 
oxygen  saturation  was  assumed  to  be  100  percent  and  the  refractive  index 
1.36  (4).  In  most  cases  the  number  of  input  photons  was  30,000. 

Optical  Parameters  Used  in  the  Model 

The  optical  parameters  such  as  scattering  and  absorption 
coefficients  S  and  K  were  calculated  according  to  Zdrojkowski  and 
Pisharoty  (41).  Path  lengths  between  scattering  events  were  based  on 
work  done  by  De  Palma  and  Gasper  (10). 

Maarek,  et  al.  reported  that  use  of  only  an  optical  collimator  is 
insufficient  to  select  only  those  photons  which  travel  near  the  axis  of 
illumination: 

It  should  be  noted  that  for  thin  samples  the  probability  of  a 
photon's  free  path  being  greater  than  the  sample  thickness  is  fairly 
high  and  so  photons  can  pass  through  without  a  single  scattering 
event.  It  is  for  these  situations  that  the  use  of  an  optical 
collimator  can  produce  good  image  resolution  using  trans illumination 
because  a  large  number  of  detected  photons  will  have  travelled 
straight  along  the  illumination  axis  (24:408). 

The  computation  of  scattering  angles  was  the  same  for  both  the 

surrounding  dilute  blood  medium  and  for  the  interior  of  the  blood  vessel 

and  was  based  on  experiments  done  by  Maarek  et  al.  on  a  dilute  blood 

sample  (24:408). 

Transmittance  and  Reflectance  Characteristics 

To  provide  a  baseline  and  normalization  factor  for  the  simulation  of 


imaging,  computations  were  carried  out  assuming  a  tissue  sample  not 
containing  a  light  absorbing  blood  vessel.  In  addition,  no  detection 
aperture  was  included  at  the  output  end  of  the  sample.  Figure  14 
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Figure  14.  Variations  of  the  number  of  transmitted  photons  (N)  for 
three  values  of  the  medium  thickness(x  );  curves  1,  2,  and  3 
correspond  to  x  values  of  40mm,  60rmi7  and  80mm,  respectively 
(24:409).  On  g?aph  (a)  the  number  of  diffusely  transmitted  photons  is 
plotted  as  a  function  of  emergence  ordinate  y  (in  mu).  On  graph  (b)  the 
same  parameter  is  plotted  as  a  function  of  emergence  time  t  (in  ps)  and 
curves  1,  2,  and  3  have  been  shifted  to  the  same  time  origin  considering 
that  time  zero  corresponds  to  the  straightforward  crossing  times  of  the 
media.  The  number  of  photons,  N,  detected  with  time  gates  of  5  ps  and 
10  ps  have  been  plotted  as  a  function  of  emergence  ordinate  y  on  graph 
(c)  and  (d),  respectively  (24:409). 
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depicts  the  results  obtained  for  tissue  thicknesses  of  40mm,  60mm,  and 
80mm.  Each  of  the  graphs  a,  b,  c,  and  d  include  results  obtained  with  no 
time  gate,  a  time  gate  of  5  ps,  and  a  time  gate  of  10  ps.  Figure  15 
represents  values  of  transmittance  and  reflectance  obtained  with  varying 
thicknesses.  Also  included  in  Figure  15  are  values  of  transmittance 
when  a  1mm  detection  aperture  only  was  used  and  when  a  10  ps  time  gate 
only  was  employed. 

Maarek,  et  al.  stated  that 

Above  sample  thickness  of  about  50mm  the  reflectance  Re  is 
approximately  constant  because  the  surface  of  the  model  has  the 
greatest  influence  on  the  intensity  of  the  reflected  photons. 
Consequently,  they  contain  no  information  about  the  deeper  regions 
within  the  model  (24:409). 

Curve  b  in  Figure  14  shows  the  flight  times  of  photons  assuming  time 
zero  corresponds  to  the  straightforward  crossing  time.  The  curve  shows 
that  as  thickness  of  tissue  increases,  an  increasing  number  of  photons 
have  flight  times  greater  than  the  straightforward  crossing  time  tg. 

In  particular,  for  a  40mm  thick  sample,  somewhat  less  than  100  photons 
of  approximately  12,000  detected  have  a  flight  time  close  to  t  .  In  the 
band,  for  a  sample  60mm  thick  or  more  virtually  all  photons  detected 
have  flight  times  greater  than  tg. 

Curve  a  in  Figure  14  indicates  that  for  a  change  from  40nm  to  80mm 
thicknesses  the  transmittance  is  reduced  from  21.4  percent  to  3.18 
percent  when  there  is  no  time  gate.  A  change  in  thickness  from  40mm  to 
80mm  with  a  time  gate  of  10  ps  causes  a  decrease  in  transmittance  from 
5.8  to  0.08  percent.  Under  the  same  change  in  sample  thickness,  a 
reduction  from  2.3  to  0.01  percent  occurs  with  a  time  gate  of  5  ps. 

The  input  power  required  to  obtain  the  transmittance  indicated  in 
Figure  14  can  be  calculated  from  a  knowledge  of  the  number  of  photons 
emerging  from  the  sample  and  the  total  number  of  photons  originally  put 


Figure  15.  Variations  of  diffuse  transmittance  T^,  collimated 
transmittance  T  .  through  a  1mm  aperture,  transmittance  T^q 
corresponding  to  a  time  gate  of  10  ps  and  reflectance  Rg  as  a  function 
of  medium  thickness  (xmax)»  The  data  are  computed  considering 
homogeneous  tissular  slabs  (24:409). 

Tp  x  10  x  =  diffuse  transmittance  through  aperture  dia.  =  1mm. 

Tjq  ■  transmittance  with  time  gate  of  10  ps. 

10  x  R  «  reflectance 


into  the  mediun.  Considering  graph  a  in  Figure  14  and  a  tissue  thickness 
of  40mm,  the  total  nunber  of  photons  received  is  the  area  under  the 


graph.  This  is 


^  -  TvuTnber  of  plv>torvs  received 

~  (1/2)  UAy 


where  N  =  400  is  the  number  of  photons  received  on  axis  and  ZAy  ~ 60rm 
is  the  width  of  the  curve  at  its  base.  Since  the  transmittance  for  this 
curve  was  21.4  percent  the  number  of  input  photons  required  is  given 
approximately  by 

r  of  inpu-t  photons 

=  (l./2)(^PO^  (GQtntnn~) 

0.214 
-  5^  lOO 


This  approximation  is  reasonable  considering  that  blood  cells  have  an 
absorption  cross  section  of  approximately  0.06  at  X=  .632^m  which 
is  much  smaller  that  the  scattering  cross  section  Os^  (  6(j^nT  at 

the  same  wavelength).  Thus,  the  probability  of  scattering  is  much 
greater  than  the  probability  of  absorption. 

From  these  considerations,  one  obtains  for  the  input  power  P: 

(£1,100  pKotyns)^.  *  to~34J-s)(4.14.\  io14HQ 

5xi-crL2,3 


5.5  5  mW 


Since  Maarek,  et  al.  assumed  He-Ne  laser  X  was  set  equal  to  .632^n.  The 


Monte  Carlo  study  also  assumed  an  instantaneous  pulse  which  was 


approximated  in  the  above  calculation  by  a  time  of  5  ps. 


Simulation  of  Imagine  by  Trans illumination 


The  following  is  a  sumary  of  the  results  which  Maarek,  et  al. 


obtained  for  two  tissue  thicknesses  and  varying  location  of  a  blood 


vessel  of  radius  R. 


The  first  case  assunes  a  40nm  thick  sample  with  a  2mm  radius  blood 


vessel.  Figure  16  depicts  the  transmittance  graphed  against  yc  which  is 


the  y  coordinate  of  the  center  of  the  blood  vessel.  The  movement  of  y 


simulated  scanning  of  the  sample  with  the  source.  All  transmit tances 


are  normalized  by  the  total  transmittance  obtained  with  no  blood  vessel 


present.  The  symbol  T^  represents  the  transmittance  without  regard  to 
either  the  time  or  position  at  which  photons  emerged  from  the  sample. 


The  symbol  Tn  indicates  the  transmittance  received  with  a  time  gate  of  n 


picoseconds  and  Tr  m  represents  a  transmittance  with  a  time  gate  of  n 


picoseconds  and  a  detecting  aperture  diameter  of  m  millimeters.  The 


detecting  aperture  allowed  for  spatial  resolution.  In  all  cases  the 


aperture  or  "window"  was  centered  on  the  illunination  axis. 


From  Figure  16  one  can  see  that  transmittance  decreases  as  the  blood 


vessel  is  moved  closer  to  the  axis  of  illumination.  Based  on  a  formula 


for  the  contrast  (C)  given  as 


•minimum  of  transmittance 
^  ”  maximum  of  iracnsiuittanc^ 


the  contrast  for  T^  was 


0.30 


1.00 


Figure  16.  Variations  of  normalized  transmittance  values  as  a 

function  of  ordinate  of  the  vessel  center  y  when  a  small  vessel 

J  c 

(R  =  2mn)  is  displaced  in  the  core  of  a  thin  tissular  slab  (X^  =  40mm 

(24:311).  Curves  are  plotted  for  values  of  transmittance  TD’  T5’  T10’ 
and  2  4  (see  text  f°r  explanation  of  symbols). 

Legend  for  graph: 
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Xc  ■  20rrm  indicates  blood  vessel  was  at  constant  depth  in  the  middle  of 
the  sample. 


With  a  time  gate  of  5  picoseconds  (without  filtering  with  an  aperture), 
the  contrast  was 


c 


5 


1  -  0.37 
1.00 


0LG3 


Similarly,  for  a  time  gate  of  10  picoseconds  (no  aperture), 


c 


10 


=  1  - 


0.52 

loo 


0.48 


The  curve  in  Figure  16  denoted  by  a  diamond  symbol  O  represents 
transmittance  values  obtained  with  a  time  gate  of  5  picoseconds  and  an 
aperture  window  of  2.4  run.  In  this  case  the  contrast  increased  to 


c 


O.Zl 

1.00 


=  0.77 


The  second  case  of  simulated  trans illumination  assumed  a  tissue 
thickness  of  80mm  and  a  blood  vessel  of  radius  R  =  4mm.  Figure  17 
depicts  the  transmittance  values  obtained  wheh  the  blood  vessel  was 
moved  laterally  across  the  middle  of  the  tissular  slab.  Contrast  values 
obtained  were  T^:  =  0.20,  T^q:  C^q  =  0.65,  T^:  =  0.50, 


0.45. 


Figure  18  represents  the  transmit tances  when  the  blood  vessel  was  at 
a  distance  away  from  the  source  of  60rrm  (20mm  from  the  unilluminated  end 
of  the  tissular  slab).  The  contrast  values  obtained  in  this  instance 
were  TQ:  Cp  -  0.11,  T1Q:  C1Q  =  0.40,  T15:  C15  =  0.37,  T2Q:  C2Q  =  0.29. 

In  addition  to  the  cases  detailed  above,  Maarek,  et  al.  also 
simulated  the  transillumination  of  an  80rm  thick  slab  with  a  2nm  radius 


blood  vessel  located  at  the  center  of  the  slab.  The  researchers 


<♦'.1 
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Figure  17.  Variations  of  normalized  transmittance  values  as  a 

function  of  ordinate  of  the  vessel  center  y  when  a  medium-sized  vessel 

J  c 

(R  =  4mm)  is  displaced  in  the  core  (X  =  40mm)  of  a  thick  tissular  slab 
(X  =  80mm)  (24:411).  c 


Legend  for  graph: 
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X^  =  40mm  indicates  blood  vessel  was  at  a  constant  depth  in  the  middle 
of  the  sample. 


Figure  18.  Same  caption  as  for  Figure  17  (24:412).  The  medium-sized 
blood  vessel  has  been  shifted  towards  the  unilluminated  side  of  the 
tissular  slab  (X  =  60inn).  Contrast  is  degraded  from  the  case  where 
Xc  =  40mm.  c 
Legend  for  graph: 


reported  that,  in  this  case,  there  was  only  a  very  slight  decrease  in 
transmittance  as  the  vessel  was  moved  laterally  with  respect  to  the 
source.  The  contrast  obtained  without  time  gating  or  spatial  filtering 
with  an  detecting  aperture  was  only  0.08.  It  was  reported  that  time 
discrimination  could  not  be  used  because  of  "uncertainties  in  the 
computed  results. "(24:412). 


Discussion  of  Results  Obtained  by  Simulated  Transillumination 

The  contrast  values  obtained  by  Maarek,  et  al.  incorporate  the  same 
general  features  one  sees  in  transillumination  of  a  breast  containing  a 
tumor  attended  by  a  proliferation  of  blood  vessels.  Notably,  Figures  16 
through  18  indicate  a  transmittance  gradient  across  the  y  ordinate  which 
to  the  eye  would  present  a  shallow  region  encompassed  by  a  light  area. 

Maarek,  et  al.  have  succeeded  in  showing  that,  for  the  optical 
parameters  used  and  the  model  chosen,  the  technique  of  time  gating 
coupled  with  spatial  filtering  by  a  detection  aperture  may  be  able  to 
increase  the  contrast  obtained  under  actual  conditions  of 
transillumination.  The  best  improvement  in  contrast  obtained  with  time 
gating  and  spatial  filtering  was  57%  over  that  obtained  without  time  and 
spatial  discrimination  (see  Figure  14). 

Maarek,  et  al.  caveat  their  study  with  the  remark  that  "...Because 
of  the  present  lack  of  experimental  data  on  the  optical  parameters  of 
biological  samples  other  than  blood,  we  have  used  a  model  with  the 
assumption  that  only  intratissular  blood  cells  interact  with  photons, 
which  could  have  resulted  in  overoptimistic  results."  (24:412). 

It  will  be  shown  below  that,  though  the  Monte  Carlo  study  gives 
somewhat  higher  transmittance  values  than  those  obtained  experimentally 
in  diaphanography ,  the  difference  is  certainly  not  so  great  as  to  throw 
fundamental  suspicion  on  the  Monte  Carlo  results.  It  will  also  be  shown 


that  the  results  obtained  theoretically  in  Chapter  II  show  a  close 
agreement  with  the  Monte  Carlo  results  described  in  this  section  when 
the  same  optical  parameters  and  thicknesses  of  tissue  are  used. 

As  an  analytical  tool,  the  Monte  Carlo  simulation  technique  has 
several  advantages  and  some  disadvantages  when  compared  to  the 
theoretical  approach  presented  in  Chapter  II.  The  simulation  allows  one 
to  keep  track  of  paths  traversed  through  the  scattering  medium  thus 
leading  to  the  possibility  of  time  domain  analysis  of  the  scattering 
process.  This  feature  led  to  the  time  discrimination  technique  which 
resulted  in  the  noted  improvements  in  contrast.  The  Monte  Carlo  method 
also  gives  one  the  capability  of  determining  the  probability  density  of 
path  times  and  lengths  from  an  output  such  as  that  pictured  on  graph  b 
in  Figure  14.  Thus,  one  again  gets  greater  insight  into  the  effect  of 
the  scattering  process  on  path  lengths.  The  theoretical  analysis 
(Chapter  II),  however,  provides  for  a  more  general  conceptual 
understanding  of  the  relationships  between  multiple  scattering  and 
transport  processes  (eg.  the  Fourier  transform  relation  between  the 
correlation  function  and  the  specific  intensity).  The  mathematical 
approach  also  allows  for  a  quicker  calculation  of  the  effect  of  changes 
in  optical  parameters,  wavelength,  and  thicknesses  of  media. 

Comparison  of  Transmit tances  Obtained  bv  Monte  Carlo  Simulation 
'  and  Theoretical  Calculation 

To  provide  some  comparison  between  the  results  obtained  by  the  Monte 
Carlo  simulation  of  Maarek,  et  al,  and  the  theoretical  approach  in 
Chapter  II,  the  calculated  value  of  the  angular  spectrum  of  the  specific 
intensity  I  (T,&  will  be  obtained  for  the  same  optical  parameters  and 
tissue  thicknesses  considered  in  the  Monte  Carlo  method. 

The  calculations  of  the  angular  spectrum  1{z0)  are  based  on  the 


formula 


I  (yfi) = I0cyp(-r)d(9)  +  I0^/(7r%r)exp  (23) 

Wo 


which  is  valid  for  large  optical  thickness T( 20: 300).  To  obtain  the 
total  intensity  received  in  27Tsteradian  in  the  forward  direction, 
l(TjS)  must  be  integrated  over  values  of  from  zero  to  71/2.  Knowing 


ds  -  sin^cos&J^ 

and  neglecting  the  small  quantity  Inexp(-7")  5(0)  one  gets  from  equation 


(23): 


I(^)  =  Wtrr/expl-Op  smzO/r]sin0ccsOS 


The  albedo  Wq  has  been  set  equal  to  unity  as  the  particles  (blood  cells) 
are  virtually  non-absorbing  at  \  =  .632^n. 

From  equation  (24), 

TX/z 

T (r:6)  -  10as/(r tt) (- l/27/Op)exp [-a^md/r] 


i-exF[-ap/T-  1] 


The  quantityO^  =  2.66(D/A)^  for  large  particles  where  D  is  the  diameter 
of  a  particle  and  X  is  the  wavelength  of  radiation.  The  scattering 


blood  cells  have  a  diameter  of  approximately  and  T  is  calculated 
from  the  formula: 

r  =  pt?z  =  H(i-h)c% 

as  Ve 

where  H  is  the  hematocrit,  Vg  is  the  volume  of  a  single  scatterer 
(blood  cell)  and  Z  is  the  tissue  thickness.  The  volume  of  a  blood  cell 
is 

Ve  =  n  (D/2  ft 

-  (7/zpmf (Zuxri) 

-  7£.97 um! 


where  t  is  the  thickness  of  a  blood  cell  which  is  approximately  7jj m. 
(19:63).  The  scattering  cross  section  CT  is  obtained  by  linear 
extrapolation  from  data  published  by  Ishimaru  (19:66)  on  the  scattering 
cross  section  of  blood  cells.  For  X-  .632^n,  0^  =  60. 

Substituting  these  parameters  into  equation  (25)  gives 

r  0.196  Io/ot 


Then,  since  intensity  is  measured  per  steradian,  one  obtains  for  a  flux 
P  in  27Tsteradians: 

?  =  0.19CI-, 


for  a  transmittance  value  of  19.6  percent.  This  value  is  in  fairly 


close  agreement  with  the  transmittance  obtained  in  the  Monte  Carlo 
simulation  of  21.4  percent.  Accomplishing  the  same  steps  for  a  tissue 
thickness  of  60mm  one  obtains  a  transmittance  of  13.5  percent  according 
to  scattering  theory.  The  Monte  Carlo  method  obtained  a  transmittance 
of  approximately  10  percent  if  it  is  assumed  that  the  number  of  input 
photons  was  30,000.  For  an  80mm  thick  sample  the  calculated 
transmittance  is  10.3  percent  while  the  reported  Monte  Carlo  value  was 
only  3.8  percent.  In  the  case  of  the  80mm  thick  sample,  the  number  of 
input  photons  must  have  been  approximately 


(1/2)  (38X60)  photons  received 
0.038 


-  30,000  joho-fcorLS 


It  is  difficult  to  assess  the  exact  nature  of  the  relationship 
between  the  theory  developed  in  Chapter  II  and  the  Monte  Carlo  method. 

Nevertheless,  it  is  clear  that  the  two  approaches  agree  fairly  well  as 
long  as  tissue  thicknesses  do  not  become  too  large. 

The  reader  should  note  that  conclusions  based  on  the  comparison 
given  here  between  the  theoretical  and  Monte  Carlo  approaches  should  not 
be  carried  too  far.  The  comparison  does,  however,  support  the  notion  of 
a  general  agreement  in  results  of  the  two  methods.  In  particular,  the 
Gaussian  approximation  to  the  phase  function  for  forward  scattering 
seems  to  be  essentially  valid  based  on  the  analysis  given  above.  Again, 
it  must  be  noted,  however,  that  for  very  thick  biological  samples  the 
results  of  the  methods  diverge. 

To  conclude  this  section,  a  comparison  of  the  transmit tances 
obtained  in  the  Monte  Carlo  tissue  model  with  those  obtained  under 

i 
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actual  breast  transillunination  would  be  useful.  Table  II  presents 
optical  densities  (38:V-9)  obtained  under  actual  transillunination. 

The  values  and  Dr  are  the  optical  densities  obtained  using 
narrow  band  filters  with  transmission  peaks  at  540  and  650  nm, 
respectively. 

The  mean  value  (X)  of  Dr  of  0.83  given  for  25  women  with  normal 

-0  83 

breast  tissue  gives  a  transmittance  value  T  =  10  *  =  14.8  percent. 

It  is  possible  to  compare  this  value  with  transmit tances  obtained  via 
Monte  Carlo  simulation  and  scattering  theory  despite  the  fact  that  the 
Monte  Carlo  tissue  model  had  only  a  0.05  value  for  the  fraction  of 
volume  taken  up  by  scatterers  (and  is  therefore,  not  a  dense  medium)  and 
yet  breast  tissue  is  relatively  dense.  Recall  that  the  equation  for 
scattering  strength 

pa  -  yr$ 

113  vc 

is  symmetric  about  its  maximum  value  of  w  =  1/2  (w  =  volume  of 
scatterers  as  a  fraction  of  total  sample  volume). 

The  equation  for/3"®  shows  that  the  values  of  w  =  0.05  and 
w  =  0.95  gives  the  same  values  of  scattering  strength.  Thus,  a  comparison 
between  transmittances  obtained  by  Monte  Carlo  simulation  and  actual 
trans illumination  of  breast  tissues  is  not  unfounded.  In  particular, 
the  approximate  agreement  of  transmittances  obtained  by  simulation, 
theory,  and  actual  trans illumination  lend  further  support  to  the 
contention  that  scattering  in  the  breast  is  primarily  forward  directed 
and  is  due  to  particles  somewhat  larger  than  red  or  near  infrared 
wavelengths. 


TABLE  II 


OPTICAL  DENSITIES 


Optical  density  values.  D 


and  Dr, 


obtained  from  diaphanograms 


of  110  women  with  breast  cancer  and  1_  women  with  benign  breast 
disorders.  The  miscellaneous  group  consisted  of  two  women  who  had  lymph 
nodes  and  one  woman  who  had  a  lipoma.  X  =  mean  value  and 
SD  =  standard  deviation. 


Histologic 

classification 


Visual  interpretation  No.  of 
of  diaphanograms  patients 


Optical  density* 


Dg  Dr 


* 

SD 

X 

SD 

Cancer 

Dark-shaded  area 

42 

2.87 

0.42 

2.03 

0.86 

Light-shaded  area 

52 

2.25 

0.49 

1.12 

0.56 

No  shaded  area 

16 

2.16 

0.58 

1.19 

0.75 

Fibroadenoma 

No  shaded  area 

37 

2.01 

0.58 

0.84 

0.53 

Mastitis 

Light-/ dark-shaded 

area  7 

2.38 

0.43 

0.93 

0.46 

Miscellaneous 
Mammary  duct 

No  shaded  area 

3 

1.69 

0.40 

0.56 

0.07 

ectasia/Flbro- 
cystlc  disease 

No  shaded  area 

82 

2.08 

0.47 

0.68 

0.24 

Normal  breast 
tissue 

*Expressed  in  i 

No  shaded  area 

density  units 

25 

2.02 

0.59 

0.83 

0.50 

(38:V-9) 
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Implementing  Monte  Carlo  Simulation  Using  the  Gaussian  Approximation 

to  the  Phase  Function. 

It  would  be  useful  to  implement  a  Monte  Carlo  simulation  of  the 
scattering  process  in  a  way  which  ties  Monte  Carlo  results  closely  to 
the  scattering  theory  presented  in  Chapter  II.  The  motivation  for  doing 
so  comes  from  the  consideration  that  the  improvement  in  contrast 
experienced  by  time  gating  and  spatial  filtering  (TGSF)  of  the 
transillumination  output  could  then  be  translated  mathematically  into  a 
reduced  optical  thickness T^.  This  reduced  optical  thickness  of  a 
tissue  medium  with  the  same  transmission  characteristics  of  a  medium  of 
thickness  where  the  scattering  cross  section  CT  is  reduced  in 
magnitude.  In  other  words,  TGSF  effectively  produces  a  "new"  medium 
with  reduced  optical  thickness  7~R.  It  is  suggested  that,  using  this 
reduced  optical  thickness,  one  could  again  apply  the  scattering  theory 
of  Chapter  II  to  determine  the  feasibility  of  mathematical  image 
processing  on  the  new  TR  medium. 

As  a  first  step  in  this  direction,  the  following  is  presented  to 
show  how  a  Monte  Carlo  simulation  can  be  created  incorporating  the 
Gaussian  approximation  to  the  phase  function  p(s)  which  describes  the 
distribution  of  scattered  radiation. 

In  capsule  form  the  steps  involved  to  implement  the  simulation  are 

(1)  Introduce  a  photon  into  the  medium  at  position  (x,y,z  =  0), 

(2)  Calculate  the  path  length  1  to  the  location  where  the  next 
scattering  or  absorption  event  takes  place, 

(3)  Determine  if  the  photon  is  absorbed  or  scattered  after 
traversing  path  length  1, 

(A)  If  the  photon  is  absorbed,  stop  and  input  a  new  photon, 

(5)  If  the  photon  is  scattered,  calculate  the  direction  (OiP)  which 


the  photon  takes,  and 


(6)  Compute  the  position  (x,  y,  z)  of  the  photon  to  track  its 
progress  through  the  medium. 

The  probability  density  functions  needed  to  calculate  1,  ,  and  at 

each  scattering  event  are 

fjd)  =  •  0<l<oo, 

Vjtp)-  )  0«P<2TT 


where  1  =  l/n(^  is  the  mean  free  path  length  (n  is  the  total  number  of 
scatterers).  Note  that  the  integration  to  arrive  at  ^6)  was  over 
limits  from  0  to Oo which  is  a  valid  approximation  since  the  contributions 

1)  and 

exponential  and  Gaussian  random  variables  would  have  to  be  generated 
from  a  random  variable  uniformly  distributed  over  the  interval  from  0 
to  1  as  generated  by  a  computer.  The  values  of  1  may  be  calculated  from 
the  equation 


to  the  integral  for  large  are  negligible.  For  p^( 


i  =  iiu[i/cu-iy] 


(26) 


where  u  is  a  uniform  random  variable  distributed  over  the  interval 
0  <  u  <  1.  The  values  of  0  can  be  calculated  from  an  approximation 
listed  by  Abramowitz  and  Stegun  (1:932). 

To  tie  the  results  of  the  Monte  Carlo  method  outlined  above  to 
scattering  theory  as  developed  in  Chapter  II,  the  output  flux  as 
predicted  by  the  angular  spectran  of  the  specific  intensity  could  be 


fitted  mathematically  to  the  output  flux  of  the  obtained  by  Monte  Carlo 
simulation  with  TGSF.  If  all  the  same  optical  parameters  are  used  in 
the  simulation  and  in  the  scattering  theory,  the  difference  in  the 
fluxes  predicted  by  simulation  and  theory  becomes  a  function  of  optical 
thickness  T. 

To  accomplish  this  mathematical  fit,  however,  the  output  fluxes  of 
the  Monte  Carlo  simulation  with  and  without  TGSF  must  be  normalized 
relative  to  each  other  to  insure  that  a  reduction  in  T ( T  ->  T  )  is 
actually  obtained.  To  see  this,  consider  typical  transmission  curves  in 
Figure  19  for  the  same  medium  with  and  without  TGSF.  If  one  does  not 
normalize  the  outputs  in  graphs  (a)  and  (b)  in  Figure  19  it  can  be  seen 
that  a  medium  corresponding  to  graph  (b)  without  TGSF  would  be  one  with 
a  larger  optical  thickness  than  in  graph  (a)  because  of  the  reduction 
in  total  transmission.  However,  one  would  expect  that  if  effective 
scattering  is  reduced  (artificially  by  TGSF)  then  graph  (b)  would 
represent  a  medium  with  reduced  T . 

Therefore,  the  areas  under  graphs  (a)  and  (b)  in  Figure  19  should  be 
normalized  relative  to  one  another  to  produce  equivalent  transmission 
values  as  shown  in  Figure  20.  Ihe  areas  under  the  graph  (c)  and  (d)  in 
Figure  20  are  equal  and,  consequently,  represent  outputs  with  equal  total 
transmission.  It  can  be  seen  that  with  transmission  normalization, 
graph  (d)  must  correspond  to  a  medium  with  reduced  optical  thickness 
because  the  spread  of  intensity  in  graph  (d)  is  decreased  with  respect 
to  that  in  graph  (c).  If  one  fixes  particle  density  f)  and  physical 
thickness  Z,  then  the  reduction  in  the  spread  of  intensity  reflected  in 
graph  (d)  corresponds  to  a  reduced  scattering  cross  sectior  produced 
by  removing  scattered  photon  noise  by  TGSF.  Consequently,  the  actual 
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Figure  19.  Typical  transmissision  curves  obtained  via  Monte  Carlo 
simulation.  Graph  (a)  shows  the  transmission  obtained  without  time 
gating  or  spatial  filtering  at  the  output  (note  the  large  spread  in 
photon  intensity^y^).  Graph  (b)  shows  transmission  obtained  with  time 

gating  and  spatial  filtering  (note  total  transmission  is  reduced  as  well 
as  the  spread  of  intensity  A y2) . 


r'.« 


it’*.*’' 


Transmission  (Nunber  of  Photons) 


— HAV2  I1 


Spread  of  Intensity  in  Output  Plane 


Figure  20.  Normalized  transmission  curves  obtained  by  Monte  Carlo 
simulation.  Areas  under  the  graphs  are  equal  (identical  total 
transmission  values).  Graph  (c)  shows  transmission  obtained  without 
time  gating  or  spatial  filtering.  Graph  (d)  shows  transmission  obtained 
with  time  gating  and  spatial  filtering. 
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optical  thickness  T 


becomes  the  reduced  optical  thickness  T^: 


(27) 


(28) 


where  CT'  <  0"s.  It  remains  to  show  how  the  value  of  CT'  may  be 


obtained  from  the  equation  for  the  angular  spectrum  of  specific 
intensity  I 


,00 


(jfi)  =  (X„/Z7r/q<^^(4»inS)e«p|-7[i-We«(p(-^]j(29) 


which  is  obtained  by  Fourier  transformation  of  the  correlation  function. 
Let  the  total  transmission  (or  flux)  obtained  without  TGSF  be  denoted  by 
,  then 

(30) 

where  N'^  (Q^ )  is  the  number  of  photons  received  in  the  output  plane  as  a 
function  of  0^  =  tan  ^(y^/z).  The  integration  is  carried  out  over  the 
angle  ,  (the  angular  extent  of  the  intensity  in  the  output  plane).  The 
subscript  1  in  equation  (30)  identifies  quantities  corresponding  to 
an  experimental  arrangement  without  TGSF.  The  symbol  y^  is  the 
coordinate  axis  in  the  output  plane  perpendicular  to  the  illumination 
axis. 

Similarly,  let  be  the  total  transmission  in  the  output  plane  with 


TGSF,  then 


where  the  subscript  2  identifies  a  quantity  obtained  with  TGSF.  If  F^ 
is  the  normalization  factor,  then 


^1  ~  FnT2(Tr)  (32) 

where  is  a  quantity  determined  by  simulation.  The  total  transmission 
(flux)  T2  has  been  integrated  over  0 ^  leaving  a  quantity  which  is  now  a 
function  of  7~^  alone.  From  equation  (29)  one  obtains 


Since  D  ,  (J  ,  and  Z  have  been  set  the  same  values  for  cases  with  and 
'  n’  a’ 

without  TGSF,  Tu  is  a  function  of  C'  (a  scattering  cross  section 

K  S 

reduced  by  TGSF).  Therefore,  equation  (33)  for  the  specific  intensity 
states  that  1(71,)  =  I((7'). 

Finally,  since  T^  =  2  1^,  and  T2  =  2  1^  CT'  may  be  determined  by 
solving  the  equation 


lz(c5')  --  Vf* 


(35) 


and  CT  has  been  set  equal  to  zero  for  a  (relatively)  non-absorbing 
medium.  Equation  (36)  can  be  solved  numerically. 

OnceO^'j  the  reduced  scattering  cross  section  is  obtained,  one  is 
then  presented  with  (or,  so  it  may  turn  out)  a  new  imaging  problem  to 
which  the  equations  for  the  limit  of  resolution  may  be  applied.  Of 
course,  the  hope  would  be  that  the  reduced  optical  thickness  7"p  is  low 
enough  that  either  a  reasonable  image  may  be  obtained,  or  that 
mathematical  imaging  processing  could  be  reapplied  with  reasonable 
expectations  of  success. 


Imaging  tumors  of  the  breast  by  transillumination  proves  to  be  a 
challenging  problem  -  and  this  for  principally  two  reasons: 

(1)  Noise  generated  by  the  scattering  process  is  not  additive,  and 

(2)  Breast  tissue  has  an  extremely  high  optical  thickness. 

Point  (1)  above  sets  fundamental  limitations  on  methodology. 

Because  the  coherent  or  average  signal  is  physically  converted  into 
noise  by  scattering,  matched  filtering  or  averaging  techniques  are  not 
likely  to  succeed  in  improving  image  quality.  A  highly  scattering 
medium  very  quickly  destroys  the  average  signal  bearing  wave  and, 
consequently,  after  a  short  propagation  distance,  there  is  little  or  no 
signal  left  to  match.  Also,  techniques  based  on  averaging  various 
outputs  suffer  from  the  lack  of  a  coherent  field  so  that  although  a  sum 
of  outputs  (perhaps  at  different  ilLumination  angles)  could  in  principal 
reduce  noise,  the  result  of  averaging  would  not  be  an  enhancement  of  the 
coherent  image  since  any  coherency  is  virtually  gone. 

Point  (2)  addresses  the  high  degree  of  scattering  which  light 
undergoes  as  it  propagates  in  breast  tissue.  The  optical  thickness 
determines  the  severity  of  the  scattering  and,  therefore,  the  minimum 
physical  thickness  through  which  one  might  obtain  a  reasonable  coherent 
image  .  For  the  aperture  and  focal  length  of  lens  chosen,  the  value  of 
for  which  a  coherent  image  could  be  expected  to  be  obtained  was  20.08 
while  for  the  biological  medium  used  in  the  Monte  Carlo  simulation  was 
nearly  750.  It  is  believed  that  the  T  for  the  human  breast  would 
likewise  be  very  large. 

Thus,  one  is  led  to  conclude  that  because  of  the  lack  of  a  coherent 
signal  in  optically  thick  tissues,  the  image  obtained  by 


trans illumination  can  be  improved  only  if  some  physical  process  to 
mitigate  the  scattering  is  first  applied.  A  technique  designed  to  do 
this  is  time  gating  and  spatial  filtering  (TGSF).  Using  TGSF  it  was 
found,  that  under  some  circumstances,  an  improvement  in  contrast  of 
77  percent  is  theoretically  possible. 

Since  the  Monte  Carlo  simulation  done  by  Maarek  et  al.  used  a  dilute 
blood  sample  as  the  scattering  medium,  it  can  be  expected  that  some 
differences  would  exist  between  the  results  obtained  by  the  Maarek  et. 
al.  Monte  Carlo  simulation  and  a  Monte  Carlo  simulation  of  breast 
transillumination.  Thus,  a  Monte  Carlo  simulation  designed  specifically 
for  breast  transillumination  could  prove  to  be  of  value. 

Such  a  study  could  be  used  in  principal  to  test  the  applicability 
various  models  and  model  parameters.  Also,  following  the  implications 
of  the  last  section  of  Chapter  III,  it  might  be  found  that  a  particular 
time  gate  and  spatial  filter  could  produce  a  reduced  optical  thickness 
Tp  for  which  mathematical  image  processing  is  a  feasibility.  Of  course, 
before  any  such  scheme  can  be  entertained  as  a  possibility,  it  is 
essential  that  experimental  data  be  obtained  for  important  optical 
parameters  of  breast  tissue  such  as  scattering  and  absorption  cross 
sections,  and  particle  densities,  concentrations  and  dimensions.  It  is 
recommended  that  if  further  research  is  conducted  on  improving  imaging 
in  breast  transillumination,  that  it  be  done  in  three  phases.  The  first 
phase  should  concentrate  on  gathering  data  as  described  above.  The 
second  phase  of  the  research  should  conduct  Monte  Carlo  simulations  to 
determine  the  reduced  optical  thickness  7"^  which  can  be  theoretically 
achieved  by  TGSF.  The  third  and  final  phase  should  be  devoted  to 
determining  the  feasibility  of  and  quantitative  efficacy  of  mathematical 
image  processing  using  the  medium  with  an  optical  thickness7^  reduced 


by  TGSF. 


Finally,  the  reader  will  recall  that  one  result  of  the  discussion  of 
the  scattering  model  in  Chapter  I  was  that  improvement  in 
transillumination  imaging  could  occur  by  simply  detecting  photons  a 
There  is  a  transmission  peak  at  this  wavelength  in  human 
collagen  of  the  dermis.  Since  the  fibrous  material  in  the  breast  is 
also  collagenous,  it  is  reasonable  to  suggest  that  experimentation  on 
breast  transillumination  at  1.2^n  would  show  increased  transmission 
through  structures  surrounding  cancers  while  yet  retaining  the  feature 
of  high  absorption  in  the  water-bearing  plasma  of  vessels  massed  about 
cancers.  An  increase  in  contrast,  then,  should  be  detectable  in  the 
image  plane.  It  is  recommended,  therefore,  that  transillumination  with 
detection  at  1.23^n  be  attempted. 
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Appendix  A:  Multiple  Scattering  Theor' 


The  following  development  of  Twersky's  theory  is  based  on  Ishimaru's 
formulation.  Consider  a  random  distribution  of  N  particles  located 
at  point  vectors  r^,  rN  in  a  volume  V.  At  this  point,  no 

specific  assumptions  are  made  concerning  the  size  or  shape  of  the 
individual  particles.  Also  consider  a  scalar  field \Jf  at  (rQ  being 
a  point  in  space  between  scatterers).  The  field  y/f  satisfies  the  wave 
equation 

(V2  +  =  O  (37) 

where  k  =  2JV  \ is  the  wave  number  describing  the  propagation 


characteristic  of  the  medium  surrounding  the  particles  (20:254). 

If ip^a  is  an  incident  wave  ("i"  meaning  incident  and  "a"  indicating 
location  rg)  in  the  absence  of  any  particles,  the  field j^/a  is  then  the 
sum  of^a  and  the  contributions  Usa  from  all  N  particles  located  at 
some  point  r  ,  s  =  1,  2, . . . ,N. 

\p  zz  (38) 

3=1 

The  symbol  U  indicates  the  field  contributions  at  a  point  r 

S  3 

(upper  index  of  U)  from  scatterers  located  at  various  points  rg  (lower 
index  of  U)  (20:254). 

The  sum  over  U  a  a)  is,  then,  the  wave  at  r  due  to  all  scatterers 

S  S  S  3 

•  s 

m  the  volume  V.  This  wave  may  be  expressed  in  terms  of  a  wa ve(p 


incident  upon  a  scatterer  located  at  rg  and  a  scattering  characteristic 
u  of  the  scatterer  located  at  r  as  observed  at  a  location  (between 

w  S 

scatterers  r  .  In  mathematical  form  U  a  is 
a  s 


a  s 


does  not 


where  uga  is  taken  to  be  an  operator  so  that  the  symbol  ug 
indicate  a  product  (in  general)  (20:254). 

The  "effective  field"  field  incident  upon  a  scatterer  located  at  a 
point  rg  is  (20:255) 

<PS  =  Ws  + 

t*l, 

t<5 

The  case  s  =  t  is  excluded  from  the  sum  in  equation  (40)  since  a 


scatterer  does  not  contribute  to  the  wave  incident  upon  it.  Figures  21 
and  22  depict  the  quantities  thus  far  described. 


1 

- > 

• 

°SS  (at  ra> 

Figure  21.  The  wave  Us  as  the  sum  of  the  contributions  from  an 

incident  wave  (p.  and  scattered  waves  from  particles  1,  2,...,  s,...,  N 
(20:255).  1 


Figure  22.  The  wave  uga(^  at  ra  due  to  t^ie  scattering  of  w ave(^f  by 
a  single  particle  S  (20:255). 


Combining  equations  (38),  (39),  and  (40)  results  in  the  following 


expression  for  the  field  : 

1  S=1  1  t-i, 


Equation  (41)  may  be  iterated  overhand  thereby  include  the  effects 


of  an  ever  increasing  number  of  scatterers.  If  such  an  iteration  is 


carried  out,  J,  becomes  (20:256) 

„  _  N  Y  N  M 


'ii^i  ♦  if t  «  (42) 

S  ws'  SJtFtW-!; 

The  first  term  in  equation  (42)  is  the  incident  wave  at  r  .  The 

3 

second  term  is  the  contribution  to^/3  from  all  single  scatterings  from 
all  the  scatterers  s  =  1,  2,...,  N.  The  third  term  represents  the 
scattering  contribution  to^f  from  all  double  scatterings.  The  terms  in 
equation  (42),  then  represent  the  scattering  contributions  from  all  the 
single,  double,  triple  and  multiple  scatterings  that  occur  in  the  random 
medium.  See  Figure  23  for  a  depiction  of  various  multiples  of 
scattering. 


(43) 


concept  of  multiple  scattering  through  one  particle  the  fourth  term  in 

equation  (42)  may  be  written  as 

_N  _N_*1  N  m  M  K  M 

Tft 

The  second  term  on  the  right  side  of  equation  (43)  has  conflated  to  a 
double  sum  since  the  sum  is  over  s  and  t  and  m  =  s  (20:256). 

Twersky's  theory  does  not  include  scattering  through  a  single 
particle  more  than  once.  As  Ishimaru  states,".  .  .Twersky's  theory 
should  give  excellent  results  when  the  backscattering  is  insignificant 
compared  to  scattering  in  other  directions"  (20:257).  Even  so,  in  other 
cases  where  backscattering  is  present  to  a  more  pronounced  degree,  the 
difference  between  the  Twersky  scattering  process  and  an  exact  process 
(including  multiple  paths  through  a  single  particle)  becomes  small  as 
the  number  N  of  particles  grows  large  (20:258). 


TABLE  III 

EXACT  THEORY  AND  TWERSKY'S  THEORY 

Shown  in  the  table  are  the  number  of  terms  included  in  the  exact 
scattering  process  and  the  Twersky  process.  As  N  (the  number  of 
particles)  grows  large,  the  difference  between  the  exact  process  and 
the  Twersky  process  becomes  very  small. 


E  =  Exact  T  =  Twersky  Difference 

E  -  T 
E 


ina 

i 

1 

1 

0 

Single  Scattering 

N 

N 

0 

Double  Scattering 

N(N  -  1) 

N(N  -  1) 

0 

Triple  Scattering 

N(N  -  l)2 

N(N  -  1)(N  -2) 

1 

N  -  1 

Quadruple  Scattering 

N(N  -  l)3 

N(N  -  1)(N  -  2)(N  -  3) 

3N  -  5 
(N  -  1) 

(20:258) 
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For  convenience,  the  mathematical  form  of  the  Twersky  process  is 


written  here  explicity. 

M  N  M 


N  N  N 


<Aa  +1  +LEl  -uNtCvf 

S--1  S -1  t=l.  S«l 


S'l 

t**  m*t, 

s 


(44) 


Integral  Formulation 

The  development  so  far  presented  is  useful  for  a  conceptual 
understanding  of  the  multiple  scattering  process  but  does  not  lend 
itself  to  use  as  a  tool  for  determining  quantities  in  practical 
problems  (20:259). 

To  make  the  theory  practical,  Twersky  derived  an  integral 
formulation  of  his  theory.  Consider  a  random  medium  consisting  of 
particles  all  of  which  have  the  same  statistical  characteristics 
described  by  a  weighting  function  w(^).  "The  variable  £  designates  all 
of  the  characteristics  of  a  scatterer  s  such  as  location  (r  ),  shape, 

o 

orientation,  and  dielectric  constant"  (20:259).  Since  all  the 
scatterers  have  the  same  statistical  characteristics,  w(^)  may  be 
written  as 


W(S)  =  w(fS)£5) 


(45) 


where ^ s  represents  all  the  characteristics  of  the  particles  with  the 
exception  of  location  (20:259). 

If  we  now  have  some  random  function  f  which  depends  on  all  the 
scatterers  we  can  express  f  as  (20:259) 


<£($)}  w •  •  *w(rsV •  -w(gvl-%’*<frrv 


(46) 


where  <f(£)>  represents  the  average  of  f  corresponding  to  scatterer 


characteristics  other  than  location.  One  can  then  write 


Nutter  of  scatterers  within 

=  VQl-uTnf  V  :dr\=  dysdysd»5 

"tot  3,1  ivuwiVer  of  st^itctcrs 

in  V 


=  (47) 

N 


wherep(rg)  is  the  nunber  of  scatterers  per  unit  volume  (20:260). 


The  average  of  f  is  (20:260) 


In  the  development  above,  it  is  assumed  that  the  scatterers  are  (to 
within  a  good  approximation)  statistically  independent  of  one  another. 
For  a  biological  medium,  investigators  assume  this  to  be  the  case  for 
wavelengths  in  the  visible  and  near  infrared.  This  assumption  may  be 
supported  by  a  consideration  of  the  far  field  of  a  particle.  For  a 
collagen  fiber  (diameters  DsslOOnm)  the  far  field  is  D  /\  =  O.O^m  for 
a  value  of  X  =  ljjn.  To  put  this  in  perspective,  a  distance  of  O.Ojr^m  is 
700  times  smaller  than  the  average  diameter  of  a  red  blood  cell.  Thus, 
the  fibers  do  not  have  to  be  a  significant  distance  from  each  other  for 
the  far  field  approximation  to  be  accurate.  In  this  case  it  may  be 
concluded  that  statistical  independence  is  a  good  approximation.  Later, 
it  will  be  shown  both  qualitatively  and  quantitatively  that  when  the 
density  of  scatterers  becomes  large  (and  the  distance  between 
scatterers  is  small),  the  character  of  the  scattering  process  changes 
significantly  (scattering,  in  fact,  decreases). 

Given  the  expression  for  <f>  in  equation  (48),  an  equation  for  the 


•i 

t'tl 


various  scattering  terms  in  Twersky  s  equation  may  be  written.  First, 
from  the  equation  for^, 

V9* (44) 

SCI  5-1  t-l,  S=l  t-l 

tils 

consider  the  second  term.  This  term  has  the  integral  form  (see  equation 


40  ) 

sV?iS>~]r 

5=1  5s 


Sfl*u.tV>i  ^dr5 


“/•WjViP^drj 


The  third  term  in  equation  (44)  becomes  (20:260) 


>"<  4‘<,  JJ  M  5  X 


Since  there  is  one  less  term  in  the  sum  over  t  as  compared  to  the  sum 
over  s,  equation  (50)  results  in 

ff  <AW>  =  ffui-rtufri  Arsi?t 


5=1  t-l, 
t*S 


r  pfcv>(rtwsd>t 


r.T.v.Ttw^ 


WT 


For  large  N  (which  is  the  case  in  breast  transillumination)  equation 
(51)  becomes 


In  general,  the  terms  in  equation  (44)  generate  multipliers  of  the 
integrals  of  the  form 

NGmXn-O  •  •  •  CM- 

M  =  - Zr*1* -  (52) 

N 

where  n  indicates  the  nth  term  in  equation  (44).  Since  the  highest 
degree  terms  in  the  numerator  and  the  denominator  of  equation  (44)  are 


lit*  M  *  1  (53) 


Therefore,  the  intebral  expression  for  is 

<iP a>  =  <P?  + 


,5dr; 


l^s  P  fP>P  C 


(54) 


Recalling  that  <1 p)~  'rtWrt  one  can  write  equation  (54) 

<i//a)=  <P±  +  |'Uj<i/'5>P<r^drs 


as 


(55) 


Iteration  of  equation  (55)  results  in  equation  (54). 

Equation  (55)  represents  the  fundamental  equation  for  the  coherent 
(average)  field  in  Twersky's  multiple  scattering  theory  and  is  the  major 


result  of  this  section. 


Twersky  s  integral  equation  for  the  correlation  of  the  field  in  a 
random  medium  at  points  and  is  consistent  with  equation  (55) 
developed  in  last  section  (20:263).  The  correlation  function  W> 
derived  by  Twersky  is  (35:99) 

<vW  <jA'  +  U\h*mz>p(^ 

where 


vs"  = 


v  will  be  referred  to  as  the  multiple  scattering  operator. 


If  equation  (1)  is  iterated  in  accordance  with  equation  (2)  the 
following  expansion  results: 

<^V*> =  </></*> + / v‘V  2p  (^dfS 

+JjjS Vs‘Vsb*',t  V»  Vm* 
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The  first  term  of  equation  (56)  is  the  pm-ik  t 
at  r  and  the  complex  conjugate  of  th<  c<>h«  r-  •  ’  • 

The  second  term  represents  the  wav*  at  r 
field  at  rg  t’.uaugh  the  mew  ha.  is- 
generated  by  the  complex  e  ■  >i  ■  j.  : *  •  •  • 
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Figure  24.  The  v  scattering  processes  giving  rise  to  the  field  at 
r  (20:264).  s 

3 

Based  on  the  development  thus  far,  the  relation  between  the 

3  3 

scattering  operators  u  and  v  becomes  clear.  In  heurisitic  terms, 

s  s 

the  scattering  operator  u  represents  the  radiation  from  r  to  r 

S  S3 

3 

through  free  space,  whereas  v  represents  the  radiation  from  r  to  r 

through  multiple  scattering  (20:270).  In  order  for  the  theory  to  become 

useful  the  mathematical  form  of  v  a  (the  multiple  scattering  operator) 

s 

must  be  obtained.  Consequently,  the  next  stage  of  the  mathematical 

formulation  of  the  scattering  problem  will  be  devoted  to  a  determination 

of  v  a  based  upon  (1)  a  reasonable  approximation  to  the  form  of  u  a 
s  s 

Q 

(the  free  space  scattering  operator  and  (2)  the  relation  between  ug  and 
vga  given  by  equation  (2)  (20:264). 


Figure  25.  Scattering  process  for  equation  (22)  (20:264) 

(a)  represents  scattering  for  the  first  term,  (b)  for  the  second,  and 
(c)  the  scattering  for  the  third  term. 


The  reader  will  recall  that  the  far  field  approximation  was  assumed 

for  the  scattered  field.  Envoking  this  same  approximation  the  wave 

..  a  „  /s,  _ _  i _ _ _ i  a _ _ /  .  nrr\ 


us  Y>  b6  wr*ttGn  as  (20:266) 


where  f(o,i)  is  the  amplitude  scattering  coefficient  for  a  single 
scatterer.  f(oYi)  represents  the  amplitude  of  the  ^component  a  field 

incident  upon  a  scatt°rer  from  the  direction  "l  and  scattered  into  the 
direction  (the  symbol  6  is  chosen  to  indicate  a  direction  of 
observation).  Given  the  expression  in  equation  (57)  and  the  relation 
between  u  a  and  v  a  [equation  (2)]  it  is  possible  to  calculate  the  form 


of  v  .  First,  from  equation  (57)  write  u  a  as  (20:266) 

O  3 


lra"  rsl 


and  by  extension,  for  the  radiation  from  another  scatterer  at  rt  to 
location  ra  one  gets 

a.  _  c  $  \ 

(s 

Figure  26  illustrates  the  physics  of  the  problem.  The  solution  to 
the  problem  of  finding  vga  consists  of  solving  the  following  equation 


Figure  26.  The  scattering  process  for  determining  the  operator  vg 
(20:269).  The  problem  consists  of  integrating  over  all  contributionssto 
vga  from  scatterers  located  at  all  positions  rt. 

To  proceed,  it  is  first  assumed  thatp(rfc)  in  equation  (57)  is 
constant.  This  assumption  can  be  seen  to  be  valid  if  one  considers  that 
over  a  given  macroscopic  region  in  the  tissue  mediun  surrounding  a 
cancer  the  density  is  relatively  isotropic  (primarily  fibrons  with 
little  fat). 

The  operator  vg^  in  equation  (2)  may  be  written  as 


-89- 


whereat  is  some  unknown  function  expressing  the  results  of  multiple 
scattering.  Then,  substituting  equation  (60)  into  equation  (2)  the 
following  equation  results  (20:269) 

d  -oo  oo 


f  ats  '  ^5  l\l(tybut  'u5>  £-tsP 


o  -oo  -oo 


The  integral  in  equation  (61)  may  be  evaluated  by  the  method  of 


stationary  phase  (20:287-291)  to  yield  (20:269) 

OO  ^  CO 
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ts 


where  r^s  |  rt  -rg|,  ^  =  jrfl  -  rt|,  andstg  is  evaluated  at  the 
stationary  phase  point.  The  contributions  from  the  regions  zt  <  zg  are 
neglected  (20:269). 

Substituting  equations  (58),  (59),  and  (62)  into  equation  (61)  one 


gets 


t  _  -  ,  Z7Tif j, 

"  1  ‘ 


(63 


It  must  be  the  case  that£as  and  ^ts  have  the  same  functional  torm  such 

t  t  I 

that 


£as  ^ts|  z 


z  .  The  following  will  show  that  the  form  for 


£ts  which  satisfies  the  condition  just  mentioned  as  well  as  equation 
(63)  is 


t  .  27Tj p£(&A)  ,, 

'  k(3-ij)  (2t  s) 


(64) 


Substituting  equation  (64)  into  equation  (63)  one  gets 

,r:,P[^as,jdt, 


=  1  +  exp 
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It  is  now  possible  to  calculate  vga  by  noting,  first,  that  by 
extension  from  equation  (60) 


(65) 


S  saS 


(66) 


Then,  substituting  equations  (58)  and  (65)  into  equation  (66)  one 


obtains 


S? 


Defining  K*  ^  '  *  equation  (67)  becomes 

K 

va  =  f (Si)  exp(iK|»&-rJ) 

5  '  HW-M 


(68) 


At  this  point  the  formulation  presented  in  this  section  has  arrived 
at  (1)  an  expression  for  the  correlation  function  (equation  (1)),  and 
(2)  an  expression  for  the  multiple  scattering  operator  v  a  (equation 
(68)). 


Appendix  C:  The  Transport  Theor 


The  transport  theory  is  described  by  a  basic  differential  equation 
which  is  equivalent  to  the  Maxwell-Boltzmann  collision  equation  used  in 
the  kinetic  theory  of  gases.  The  usefulness  of  the  theory  is  not 
limited  to  photon  transport  but  has  been  applied  to  many  physical 
phenomena  ranging  from  underwater  visibility  and  neutron  transport,  to 
the  propagation  of  radiant  energy  in  the  atmospheres  of  planets,  stars, 
and  galaxies  (19:148). 

The  fundamental  quantity  in  the  transport  theory  is  the  specific 

intensity  I(r,'s)  (also  called  the  radiance).  I(r,'s)  is  measured  in 
-2  -1  -1 

Wm  sr  Hz  and  is  the  average  power  flux  density  within  a  unit 
frequency  band  centered  at  frequency  V  and  radiating  in  the  direction  s' 
from  a  point  r.  See  Figure  27. 


The  amount  of  power  dP  flowing  within  a  solid  angle  dCJthrough  an 
elementary  area  da  (with  vector  s'q)  in  a  frequency  interval 
V  to  d^  is  (19:149) 


dP=  Kfs)  cos0d  1  d6J  d  V  (w  a±t  s) 


(69) 


It  is  also  convenient  to  define  the  average  intensity  U(r)  in  terms 


of  the  specific  intensity: 

U(?  )  =  (!/%/ 

471 


where  the  integration  is  taken  over  all  47Tof  solid  angle  (19:152). 

If  the  specific  intensity  is  independent  of  s',  the  radiation  is 
isotropic.  On  the  other  hand,  if  I(r,s)  depends  on  $,  the  specific 
intensity  scattered  bypds  particles  in  a  volume  ds  can  be  expressed  in 
terms  of  the  amplitude  scattering  coefficient  by  (19:156) 


Figure  28.  Scattering  of  specific  intensity  incident  from  direction 
s'  into  direction  ^  (19:156). 

where  fCs/s')  indicates  the  amplitude  of  radiation  incident  in  direction 
"s'  and  scattered  into  direction  's. 

It  is  convenient  to  define  the  phase  function  p(s,£')  according  to 


the  following  equation: 


This  phase  function  provides  a  quantitative  description  of  the 
directional  scattering  characteristic  of  a  single  scatterer  (from  s' 
into  s ).  It  represents  the  amount  of  scattered  power  in  a  particular 
direction  and  has  no  relation  to  the  phase  ol  the  wave  (the  terminology 
"phase  function"  has  its  origin  in  astronomy).  The  albedo  of  a  particle 


0’s 

wn=  --- 
0  rr 


(where  O’  =  O’  +  O’ )  is  dependent  on  the  phase  function 

U  S  a 


(19:11): 


Wo  =  <V«r)f  P^s,s')dw  (72) 

1  -47T 

Using  the  terminology  and  symbology  just  presented  the  differential 
equation  for  the  specific  intensity  takes  the  form 

-  -PC!-#.*)*-  f  (73) 

\ 7T 

in  the  absence  of  sources  inside  the  scattering  medium  (19:157).  In 
general  a  source  term  €(r,s)  would  have  to  be  added  to  the  right  side  of 
equation  (73). 

Consider  the  terms  in  equation  (73).  The  first  term  on  the  right 
side  of  the  equation  is  due  simply  to  the  losses  in  a  volume  ds  from 
scattered  and  absorbed  power  and  has  the  form  of  the  Lambert-Beer  law. 
The  second  term  represents  the  added  contributions  to  the  specific 
intensity  from  all  scatterers  in  the  volume  ds. 

It  is  possible  to  write  equation  (73)  in  terms  of  a  nondimensional 


"optical  thickness"  defined  by 


whereupon  equation  (39)  becomes  (19:157) 


=  -l(r,t)+(l4r)|  P^.^)l(r,^75) 

The  optical  thickness  will  be  used  extensively  in  sections  to  follow 
and  in  the  calculation  of  the  limit  of  resolution  of  an  imaging  system 
in  the  presence  of  a  scattering  mediun. 

The  integral  form  of  equation  (73)  will  be  used  in  the  following 
section  to  derive  the  relationship  between  the  analytic  and  transport 
theories.  Therefore,  the  integral  equation  based  on  the  specific 
intensity  is  presented  here.  Its  derivation  is  provided  by  Ishimaru 
(19:161-162).  After  integration  over  the  total  volume  V  of  the  medium, 
the  following  equation  is  obtained  (19:162): 


U(fV>  =  U„i(raH  fcU&ir)  f?(i,Sr)l(r%t)dLO'] 


An 


*  exp[-|ra-f||] 
47T|f|-r5|2- 


(76) 


where 


U(fa)  =  average  iTvtensity 


=  (4471-)/  I(*V,$>dcJ 

V 


(77) 


and 


a  •  (r^  =  average  reduced  intensity  =(%zr) 

ri  r  J  ^47T 

=  (  Mtt)  / fa. ,  s  )  ex  P  (- 1  r^-  r0| )  du 
t7T 


(78) 


where  I  ^  is  the  reduced  intensity  which  is  the  result  of  power  losses 
from  the  incident  intensity  I^(r^,'s)  (r^  =  point  of  incidence  on  the 

medium)  due  to  absorption  and  scattering  in  accordance  with  the 
Lambert-Beer  law  of  propagation.  This  reduced  intensity  propagates 
along  a  straight  line  connecting  Tq  and  tq,  the  point  of  interest. 
Figure  29  illustrates  the  quantities  involved  in  equation  (76). 


Figure  29.  Physical  quantities  appearing  in  equation  (76)  (19:162). 


The  first  term  in  equation  (76)  is  the  average  intensity  which  results 
from  the  intensity  reduced  in  strength  and  by  losses  due  to  scattering 
and  absorption.  The  second  term  represents  additive  contributions  to 
the  average  intensity  at  r  due  to  radiation  scattered  to  that  point 
from  all  scatterers  in  the  volume  V.  Note  that  these  contributions  due 
to  scattering  are  each  diminished  by  spherical  expansion  (indicated  by 

..ue  (Vi  factor)  and  by  attenuation  (indicated  by  the 

I  r»l 

ex  Pl+a  -  rsl  ]  factor)  (19:162). 


From  the  correlation  function  of  analytic  theory  the  intensity  may 


be  written  as 


<|<aT>  =  +  /isY^s|V(^drs  (79) 

Vsa=-uf  +  [u^vsP(f%)d^t  (80) 
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where^/  '  has  been  set  equal  to^/  to  calculate  the  intensity  at  rg. 
Again,  the  form  of  vg  is  given  by 


expCiK  lra-rs  1) 
r-t.  -  rsl 


(68) 


where  K  =  k-f  [_ ($J o)/jc J  and  l  is  a  unit  vector  in  the  direction  of 
the  incident  wave.  The  term  j  vsa|2<|^f  p>  i-n  equation  (79)  represents 
the  contribution  to  the  intensity  from  all  incident  directions  i. 

To  demonstrate  that  the  correlation  function  and  the  specific 
intensity  constitute  a  Fourier  transform  pair,  it  will  be  assumed  that 
the  Fourier  relation  exists  and  that  the  equation  for  the  average 
intensity  (equation  (76))  follows  from  this  assumption.  Consequently, 
the  correlation  function  takes  the  form  (20:275) 

<l//(4)l//Vo>=  r(ilA>  r(r,r^=  fl^,S)exp(iM4)dw 

where  ?=  i/z  (4  6  and  Q  =  YL-rb  ,  and  Kf  is  the  real  part  of  K 

given  in  equation  (68). 


The  intensity  at  a  point  r  can  now  be  written  as 
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r: 
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<|i/'Ct)|z)-  T(r,o)  -  Jl(r,  s)dw 


Following  equation  (81)  and  operating  on  <|^/S|^>  with  | vga [ ^  the 
equation  for  the  scattered  wave  from  a  particle  may  be  derived.  First, 
the  operator  }v  a|2  is  calculated: 


|va'J2  _  V  V“*  -  I  f(s  s')l2  exp(-iK  lirv-f^t) 

5  1  l4-r5|  l^-r5| 


Using  the  relation  2  ImK  =  i'1  (K-K*) r  pcrt  (Im  denotes  "imaginary  part 


of")  equation  (82)  becomes 


[£(£sf)|2  P (  P°t\  ra~ rs  1 ) 

*  lraL-r5l 


where  o  and  i  have  been  replaced  by  "s'  and  s'* .  The  equation  2  ImK=p  CTt  is 
a  result  of  the  forward  scattering  theorem  (19:14).  Having  obtained 
vgaf  the  quantity  (jvsaj2,  may  be  evaluated  to 


„V/l,/,s|z\-LW3$'Mz  rs|)  AfS 

vsl  I  /-f^5'5  A  — 


Now  using  the  equation  for  the  phase  function  pCs,^ ): 


P(s, S')  »(4ff/^)lf^|a- 


and  the  average  intensity 


(86) 


U(ri')  -  (V47r)f  Itf$)c!u> 
-4 TT 

the  integral  in  equation  (79)  may  be  written  as 


As  for  the  first  term  in  equation  (79)  it  should  be  recognized  that 
the  coherent  intensity  ¥\2  >  attenuates  in  just  the  same  manner  as  the 
reduced  intensity  incident  on  the  scattering  medium.  Thus,  the 
following  equation  is  finally  obtained 


which  is  identical  to  the  equation  for  average  intensity  arrived  at  by 
the  transport  theory.  Equation  (87)  was  obtained  by  making  the 
assumption  expressed  by  equation  (3)  which  relates  to  the  correlation 
function  to  the  specific  intensity.  Thus,  equation  (87)  establishes 
the  result  that  the  correlation  function  and  the  specific  intensity 
are  related  by  a  Fourier  transform. 


Appendix  E:  Plane  Wave  Illumination  of  a  Plane-Parallel  Medium 


A  critical  question  to  ask  is  to  what  extent,  if  any,  it  is  possible 
to  use  image  processing  techniques  to  obtain  or  improve  an  image  of  a 
breast  cancer.  As  an  approach  to  this  question,  the  problem  of  plane 
wave  illumination  of  a  mediun  bounded  by  parallel  planes  oriented 
perpendicular  to  the  direction  of  incident  radiation  will  be  considered. 

The  first  step  towards  a  solution  to  the  problem  of  plane  wave 
incidence  is  to  calculate  the  specific  intensity  using  the  transport 
theory  under  what  will  be  referred  to  as  the  small  angle  approximation. 
In  this  approximation,  a  single  scatterer  is  considered  comparable  to  or 
larger  than  the  wavelength  of  illuminating  radiation  and  so  scatters 
incident  radiation  in  primarily  the  forward  direction.  This 
approximation  is  in  accordance,  then,  with  the  model  presented  in 
Chapter  1. 

Once  the  specific  intensity  has  been  calculated,  the  result  of  the 
last  section  will  be  envoked  in  order  to  obtain  the  correlation  function 
and  finally,  the  limit  of  image  resolution. 

Using  the  relation 

S.gradKr.S)  (88) 

a  S 


the  equation  of  transfer  (transport  equation)  may  be  written  as 

S.gttJ  ItfS)  r  ^  f (89 

47r47T 

in  the  absence  of  a  source  inside  the  medium  (19:235). 

In  Cartesian  coordinates  r  *  xx  +  yp  +  zz  where  x,y,  and  “'z  are  unit 
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vectors  in  the  x,y,  and  z  directions,  respectively.  The  direction  s 
(the  direction  of  scattered  light)  in  terms  of  direction  cosines  (l,m,n) 
is  given  by  (19:234) 


where 


A  i  a  A  A 

S  -  1*  +  TY\y  -h  Ti  z 

1  -  5in0cosV? 
mn  -  sinSsin^ 


n=  cos  6  (93) 

Also,  the  differential  element  dldm  which  is  set  equal  to  ds  where 
S  =  K  +  -my  has  the  form 

ds=  dldm  =  cosfisin^dSdV5  -  ndcj  (94) 

Under  the  small  angle  approximation,  the  angle  should  always  be 

small  and,  therefore,  n  =  cos  1.  Further,  although  the  limits  of 

2  2 

integration  on  1  and  m  should  be  1  +  m  1,  the  contribution  to  the 

2  2 

specific  intensity  from  regions  where  1  +  m  1  is  small  and, 

therefore,  the  limits  of  integration  on  1  and  m  can  be  extended  to  ±00  : 

f doj  fd  1  fdrt\  =  /ds  (95) 

-W  loo  loo  1 

It  is  also  assumed  that  the  phase  function  pCs,s')  is  a  function  of 

/\  A  f 

s  -  s  . 

Using  the  above  approximations  the  transport  equation  becomes 


V » '  • 
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where 


f  -  +  i  i  Vt 


-  HI  +di  *.  K  a 
'  dx  dy  >  s=Ut"y 


The  quantityPn  designates  the  number  density  of  scatterers. 

At  z  =  0  (the  boundary  of  the  medium)  the  specific  intensity  is 

I(OAS)  =  (97) 

The  general  solution  to  equation  (96)  was  obtained  by  Ishimaru  using 
Fourier  transform  method.  The  detailed  solution  to  this  problem  is 
presented  in  Appendix  J.  The  final  solution  for  the  specific  intensity 
is  given  by  (19:238) 


I(?  -is‘q)F0^,5+-A:z)K(zA,5) 


(98) 


where 


K(2,k,  f)  =  exp[-|/%r|  i-  PC?4*  (z  -2')jjdzr] 

where  q  -  (P^~PX)  - 


(99) 


The  vector  K  is  a  spatial  frequency  and  the  Fourier  transform  of  the 
specific  intensity  with  respect  to  p  =  x£  +  yy  is 

■/{ife,/3,  ?,>}  *  s  )= /" [  (2, A  s)exp(i<-p)clp  (100) 

where  7H  indicates  the  Fourier  transform  operation.  The  inverse 
Fourier  transform  of  q(z,K,s)  is 


I(zA§)=  ^”2. j !(*,*,  s)exp(i**^)<i* 


(101) 


In  addition,  the  Fourier  transforms  of  the  phase  function  p(s)  and  the 
flux  density  are  given  by 
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P(c|)  = 
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and 
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(103) 
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(for  the  definition  of  ^(z,  Kt  see  Appendix  J). 

Consider  now  the  case  of  plane  wave  incidence  on  a  tissue  medium 
bounded  by  parallel  planes  located  at  z  =  0  and  z  =  d.  At  the  boundary 
z  =  0  the  incident  specific  intensity  can  be  written  as 


l0(As)  =  Ic,6(S) 


(104) 


and  the  flux  at  z  *  0  is  (see  Appendix  J). 

F0(K,q)=  (znf  108(K) 


(105) 


To  approximate  forward  scattering  a  Gaussian  form  for  the  phase 
function  may  be  assumed  whereupon  the  phase  function  p(s)  becomes 


-p(s)=  ^V^expC-OpS2-) 


(106) 


2 

where  Q^is  proportional  to  (D/^;  and  D  is  the  particle  diameter  and 
is  the  wavelength  of  illuminating  radiation. 

If  equation  (105)  and  the  Fourier  transform  of  equation  (106)  are 
substituted  into  equations  (98)  and  (99)  respectively,  one  gets  for 
K(z,A(,q)  and  I(z,^,s)  the  following  solutions: 


K(z,k,3.)=  exp 
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and 
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The  steps  that  lead  to  equations  (107)  and  (108)  may  be  found  in 


Appendices  J  and  K.  Appendix  J  presents  the  general  solution  of 


equation  (96)  leading  to  equation  (108).  Appendix  K  describes  the  steps 


required  to  calculate  the  Fourier  transform  of  the  phase  function  to 


arrive  at  equation  (107). 
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Appendix  F:  The  Correlation  Function  in  the  Small  Angle  Approximation: 

Plane  Wave  Solution 


The  general  solution  for  the  specific  intensity  is: 

!(*,  ,?)  [ d*d?e*P(- i*P+is^F0(^«)R(e/,4) 


(98) 


where 


f0  (*,<*)  =  (ho0,  ^)exp(-i^/3+is-^')dpds 


K(^,^4')  =  ex 


P($  =J?C<)ex^ds 

For  plane  wave  incidence,  the  specific  intensity  at  t  *  0  is 

I0(M  =  I05(s) 

and,  therefore  the  flux  density  at  z  =  0  becomes 

Czr>%6(K) 


(109) 
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To  approximate  forward  scattering,  the  phase  function  is  written  in 
Gaussian  form: 


■pcs')  =  4o^/0exp  (- 

where  is  proportional  to  (D/X)^  and  Wn  *0"s/^t.  By  equation  (111) 


(A) 


P(q)  becomes 


P(q}=  47rW0exp 


(in) 


The  Fourier  transform  calculation  to  derive  equation  (111)  from 
equation  (4)  may  be  found  in  Appendix  K. 

Substituting  equations  (112),  (113),  and  (111)  into  equation  (98)  one 

gets 
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-  ex  p ( 1/c-p )  K  (z,  q)  ^  els  ‘  ^ ^ 


(115) 


Evaluating  equation  (115)  for P=  0  yields 


^JK(z4)CiH^  (H6) 

Since  the  correlation  function  and  the  specific  intensity  are  Fourier 
transform  pairs  and  equation  (116)  is  an  inverse  Fourier  transform,  it 
follows  that  f  * 

=  ^{u  <*11) 

-  I0R(Z,9^  (117) 


However,  by  equation  (106), 


K («,€')=  exp[-J^[l-  4fP(?)]dz'l 

=  expf-/^rtz.[l-Wcexp(- 

Therefore,  by  equation  (117)  the  correlation  function  becomes 

rO,q)  =  I„e  X  p  [-P  cr  2  [l  -  wcex  p(-  qV^)]  ] 


* 


where 


r  (zfij) = rcE,^| 


Appendix  G:  Image  Resolution  Limit  in  a  Random  Medium 


Consider  a  monochromatic  plane  wave  normally  incident  upon  a  medium 
consisting  of  randomly  distributed  scatterers  with  absorption  cross 
section  Oa  and  scattering  cross  section  and  optical  thickness 
+  CT  )*d  where  d  is  the  physical  thickness  of  the  medium. 


Figure  30.  A  monochromatic  plane  wave  propagating  through  a 
random  distribution  of  scatterers  (20:301).  The  resultant  scattered 
wave  is  focused  onto  the  x,y  plane  (Fourier  transform  plane).  The  lens 
has  radius  a  and  focal  length  /. 

The  resultant  scattered  wave  is  focused  by  a  lens  of  radius  a  and 
focal  length  f  onto  the  x,y  plane  as  shown  in  Figure  30.  Since,  in  the 
absence  of  scattering,  the  lens  would  produce  the  Airy  pattern  in  the 
x,y  plane  (representing  the  diffraction  limit),  the  actual  intensity 
distribution  that  results  from  scattering  can  be  compared  with  the 
iffraction  limited  case  to  get  a  sense  of  attainable  resolution  for  an 
arbitrary  object. 

In  the  x',y'  plane  (the  plane  where  the  lens  is  located)  the  field 
is  \p(z,p')  ip'  -  x'x  +  y'y).  The  field  at  the  focal  plane  (x,y)  can  be 
written  in  general  terms  using  the  Kirchoff  formula  (20:302): 


%{p) = 4  / 

S 


(120) 


where  the  integration  is  taken  of  the  area  S  of  the  lens  and  r  is  the 
distance  between  the  point  on  the  aperture  (lens)  and  the  point  on  the 
focal  plane  (2:302). 

Using  the  Fresnel  approximations 


^  ^  (l/f  1  explk[f  V  cx'-xfH/-yf  ] 


and 


~  -  —  (x'2+y'2) 


li 


and  substituting  these  approximations  into  equation  (116)  one  obtains 


1 P((P)  =(k/2TTfi)ext.(ikf  + 

JS 


(121) 


where J$  =  xx  +  yy  and  0  =  x'x  +  y'^. 

To  get  the  intensity  in  the  focal  plane  the  quantity 
PM1  must  be  calculated.  This  is  done  in  the  following  manner: 


P{(P)PtM  ~  (ikf  + 

5 


X  {Ml  )*e  *?  c1  -  jfM  6  jO()P&,p’)  <k 

Z  Z  (122) 


Carrying  out  the  indicated  multiplication  results  in 


(123) 


!<405)|z=  'Jfextf.  f 

3  _  ,  ^ 

If,  following  Ishimaru  (2:302),  the  identifications/Od'  =P^'  ~P 2’> 
4'  =  1/2  (/^  +/^2'),  and  dp'd<^'c  =  c^cp2  are  made,  equation  (123) 
may  be  written  as 


(124) 


Note  that  the  equation  is  not,  strictly  speaking, 

mathematically  true.  However,  under  the  integration  indicated,  the  end 
result  is  the  same  whether  one  uses  A' $2  or  A'A'  since  all 
points  in  the  aperture  are  included  in  either  case. 

Using  the  result  obtained  in  the  last  section  for  the  correlation 
function,  namely  that 

ra,/3d')=  <<A(z,4’)^,4')> 

=  Ioext.{-r[i.w0exf(-^)]|  (125) 

i? 

where  7"  =  PnC^z,  it  can  ^  seen  that  equation  (124)  takes  the  form 

P{09>  =  <|#5lZ> 


f\ jN  <r  %£ 

J 


(126) 


The  quantity  P^(p)  is  called  the  point  spread  function,  and  in  the  case 
of  plane  wave  illumination,  is  the  intensity  in  the  focal  plane  as  the 
result  of  imaging  a  point  source  at  a  large  distance  from  the  scattering 


medium.  An  evaluation  of  the  point  spread  function  provides  information 
as  to  how  the  scattering  (random)  medium  "breaks  up"  the  forward 
propagating  plane  wave  into  a  myriad  of  phase  fronts  propagating  in 
various  and  constantly  changing  directions  as  the  wave  strikes  particle 
after  particle  in  the  medium.  It  can  be  seen  from  equation  (126)  that 
the  correlation  function  T  determines  the  behavior  of  the  point  spread 
function  and  is,  therefore,  from  a  mathematical  point  of  view,  the 
quantity  representative  of  the  scattering  process. 

To  carry  out  the  integration  indicated  in  equation  (126)  it  is 
useful  to  define  a  function  F (0)  (20:302): 


Y0)  = 


?  Ip  I  <  a. 


(127) 


5  IP  |  >  SL 


Given  equation  (122),  Pf(p)  becomes 


PfP  W 

~Oo 

Using  the  relation  Pd'  =P^'  ~  an<^  ^c'  =  1/2(0^'  +  @2"* ’  t*ie 

-J* 

integration  of  WtPj' )  can  be  written  as 


(128) 


“(Pd) = jfe(Pc,+^(Pc  -  ^ 


(129) 


The  quantityK  (P^t )  is  the  area  of  the  intersection  of  two  circles  of 
radius  a  and  whose  centers  are  separated  by p^i .  The  solution  to 
equation  (129)  can  be  obtained  geometrically  resulting  in  (20:303) 
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K(Pd)  =  K(^V- 


>2a  (i30) 


Now  equation  (128)  may  be  written  as 


?f(P)  =Pf(p) 


~<ftf  f  “teWrvimVi 


(131) 


Note  that 


AP,'  =  (x$+y9)*  [0<i-xz)$  <- (yi-y^y] 

=  X(*i-xi'Hy(y1-yi) 


Transforming  equation  (129)  to  polar  coordinates  by  the  relations 

x-pc  os<p  y=  Psin(p 

*1=  P£°3@  yx=  Qsi-nO 

Xz=  PZC0$6  yt-  pcos9 


(132) 


dp*)  ip[  =  fi'^0 


(133) 


and  substituting  equation  (133)  into  equation  (132)  one  obtains  after 


algebraic  manipulation: 


P'P^  -pPd[cos(#-<p)] 


(134) 


Then  equation  (131)  becomes 


Pf(  ]  ~(z  ifj 


P'd<ip’A$  e*f>£f^cos(0-<$r(2,  pJ)K  CiOjO 


(135) 


Upon  interchanging  the  order  of  integration  one  gets 

2  2&f  27T 

^(hf  exp[f 

^  ro 

And  using  the  identity 


J0(a)^  ^/e^[-jacosC^-<^)]d0 


Pf^0)  can  be  written  as 


?tW'£tr  J  c^pd'Jo(|-^’)]r(z^)K^) 


and,  therefore, 


p^rcz^M(pp 

JQ 


where  is  given  in  equation  (125)  and  Kjpd’)  is  given  in  equation 

(130). 

Equation  (7)  (the  point  spread  function)  gives  the  intensity  in 
the  focal  plane.  Note  that,  in  the  case  of  no  scattering, 

A,  in  =yOpj.z  zero  and  I” ( z ,/^j »  )  (equation  (125))  becomes  equal 
to  Iq.  Therefore,  in  the  absence  of  scattering,  from  equation  (126)  the 
intensity  in  the  focal  plane  is  (20:302) 


?aCP)  a  JSf>  ^pp’tf&p’i 


(136) 


k  f 

Letting  r'  =  -jyj  one  obtains  dn‘  =  — dr*  and  equation  (136)  becomes 

^  z  (K/£)pa k 

PaC  }  Jo  Ci"r)£  <*r'f 

&7T-£)2L  J  k  0  J 

r(V£)Pa 

=  J0Cr')Jr' 


Then  using  the  relation 

X 


£  Ic£)  df  =  x  Jjfx) 


one  obtains 


?a(  )  = 


-  J„£: 


=  i0e^)jt2(^) 


(137) 


Equation  (137)  is  the  Airy  pattern  created  by  a  lens  of  radius  a  and 
focal  length  f. 

Equation  (7)  is  the  major  result  of  this  section  and  will 
constitute  the  basis  for  the  following  analysis  of  the  limit  of 
resolution  one  can  expect  from  an  imaging  system  in  the  presence  of  a 
scattering  medium. 

To  arrive  at  a  useful  interpretation  of  equation  (7),  it  is 
convenient  to  consider  the  case  of  large  optical  thickness  (T  >>  1) 
corresponding  to  a  situation  in  which  there  is  strong  scattering  (large 
^s)  and/or  a  large  physical  thickness.  In  this  case,  the  correlation 
function  I”*  takes  the  form 


~  ('r>  l>-  w°+  w°(q2/4at))] 


(138) 


where  the  approximation 


exf(-az/4)~  1-  (<f4) 

has  been  used. 

To  equation  (138)  may  be  added  lQexp(-7")  since,  if  T  is  large,  this 
additional  term  can  contribute  only  a  small  amount  to  the  value  for  I"!, 
If  these  approximations  are  made,  equation  (125)  becomes 


r  I0exp(-T)  |('r)  [i-vv0+w0 


(139) 


With  r expressed  in  this  form,  it  becomes  clear  that  the  first  term 
(iQexp(-T))  represents  a  coherent  component  of  the  correlation  function 
propagating  through  a  medium  in  accordance  with  the  Lambert-Beer  law. 

The  second  term,  then,  represents  an  incoherent  component  of Tindicating 
the  presence  of  scattering.  By  equation  (139),  then,  it  has  become 
possible  to  seperate  Tinto  coherent  and  incoherent  parts  (20:304). 

Given  the  form  of  Tin  equation  (139)  and  considering  equation  (137) 
(the  Airy  pattern  under  no  scattering)  it  can  be  seen  that  equation 
(7)  has  a  solution  for  large  T  which  is  the  sum  of  coherent  and 
incoherent  components. 

Accordingly,  replacing  Iq  with  iQexp(-T)  in  equation  (137)  one  gets 
the  coherently  propagating  portion  of  P^(p): 

P c<p)  =  coherent  component  of  the  intensity 


-  Ioexp(-T)  (SL/p')2  Jl 


(9) 


The  incoherent  part  of  the  intensity  is  unchanged  from  its  form  in 


equation  (7)  (in  accordance  with  the  approximation  indicated  in 
equation  (139),  and  therefore 


P^(p)=  incoherent  component  of  the  intensity 


~2 


(140) 


where 


=  wf-' r[i-  w0+we 

1  4a*  1 


(141) 


and  q'  =  Kr/l^i .  Using  equations  (9)  and  (140),  the  intensity  in  the 
focal  plane  can  be  written  as 


Pf(p)  =  Pt(P)  +1j(P) 


(8) 


At  this  point  it  is  convenient  to  define  the  correlation  distance-  a 
concept  which  will  be  employed  later  to  recast  equation  (140)  into  a 
form  which  lends  itself  to  gaining  a  clearer  understanding  of 
limitations  on  resolution. 

The  correlation  distance  Pq  is  the  distance  from  the  axis  of 
illumination  at  which  the  correlation  function  fall  to  e  ^  of  its  value 

at  Pd  =  0. 

The  incoherent  portion  of  the  correlation  function  I"?  (equation  139) 

at  Pd  =  0  is 

Ii<V2|)=  I0exp{-T(i-W0)J  ««> 

Taking  e-*  of  I"i  one  obtains 


a  x 
KtPo. 


[  T0  exp  (-TCi-Wo)]  e  1  =  I0e*f?E (143) 


where  has  been  set  equal  to Pq.  After  eliminating  the  factor 
lQexp[-T(l  -  Wq)]  from  both  sides  of  equation  (143)  one  arrives  at  the 
results  (20:304) 


1/2 


(144) 


It  may  be  assumed  that,  due  to  intense  scattering  (large 7),  the 

correlation  distance  in  the  case  of  transillumination  is  much  smaller 

that  the  size  of  the  aperture.  Under  these  circunstances  ,  K(Pd')  in 

2 

equation  (140)  may  be  set  equal  to  K(0)  =7Ta  and  equation  (140)  may  be 
integrated  from  0  to®.  Consequently,  (20:304) 

?i(  )  =ik1/Pi^dJ°(^Pi)ri(Z^a)7ra2  (145) 

This  equation  may  be  solved  using  the  integral  formula 

V 

/V+l  J.  o  2 

J^(at)exp(-fV)t  dt  =  - exp  (-^ 

to  get 

Pi  CP)  =  a^*>)[VrW°Kr 

Then  defining  the  quantity  (20:304) 

t 

one  obtains 


?i(p)  =  ioWPif^ptni-^-ip/Pif] 


do) 


Using  equations  (9)  and  (10)  it  is  possible  to  gain  a  clearer 
understanding  of  the  imaging  process  in  the  presence  of  scatterers. 
Equation  (9)  indicates  that  the  coherent  portion  of  the  intensity 
propagates  with  the  same  pattern  as  in  free  space  except  for  the 
Lambert-Beer  type  of  attenuation.  According  to  equation  (10),  on  the 
other  hand,  the  incoherent  intensity  spreads  out  in  the  focal  plane. 


Appendix  H:  Phase  Retrieval 


The  phase  retrieval  problem  is,  in  general  terms,  concerned  with  the 
question  of  how  spectral  information  may  be  derived  from  a  knowledge  of 
only  spectral  magnitude  (or  intensity)  information.  The  importance  of 
phase  in  imaging  is  well  attested.  From  a  complete  knowledge  of  the 
Fourier  transform  (including  both  magnitude  and  phase)  it  is  possible  to 
determine  by  inverse  Fourier  transformation  the  image  of  an  object. 

This  fundamental  capability  is  guaranteed  by  the  uniqueness  of  the 
relation  that  exists  between  a  function  f(x'*  and  its  Fourier  transform 
p(£).  However,  if  it  is  the  case  that  only  the  magnitude  of  the 
Fourier  transform  is  known,  it  is  not  possible,  in  general,  to  infer 
uniquely  the  characteristics  of  an  object  f(x)  as  any  number  of  objects 
can  have  the  same  Fourier  magnitude.  The  distinguishing  feature  among 
these  objects  in  the  spectral  domain  lies  in  the  phase  of  the  Fourier 
transform. 

In  recent  years  some  important  advances  have  been  made  regarding  the 
problem  of  reconstructing  an  image  from  only  the  knowledge  of  the 
Fourier  magnitude.  These  methods  employ  phase  retrieval  algorithms  when 
theoretical  foundation  is  based  on  the  fact  that  the  real  and  imaginary 
parts  of  an  analytic  signal  are  Hilbert  transform  pairs  (16) .However,  it 
is  not  the  purpose  of  this  research  to  explore  the  theory  of  phase 
retrieval  as  such.  The  interest  here  is  to  employ  the  scattering 
theory  of  Chapter  II  to  analyze  the  efficacy  of  phase  retrieval  as  it 
might  be  applied  to  the  problem  of  image  reconstruction  in 
t rans i 1 lunina  t ion . 

For  the  purposes  of  analysis,  the  phase  retrieval  problem  to  be 
considered  is  characterized  as  the  problem  of  reconstructing  an  image 


from  a  single  intensity  measurement  in  the  Fourier  transform 
(diffraction)  plane  and  an  estimate  (based  on  theory)  of  the  point 
spread  function  of  an  imaging  system  in  the  presence  of  a  scattering 
mediun.  The  particuclar  phase  retrieval  algorithm  that  would  be 
exercised  on  this  is  assumed  to  be  of  the  error-reduction  type  such  as 
the  Gerchberg-Saxton  algorithm.  The  error  reduction  method  was  chosen 
for  analysis  because  it  has  the  property  that  the  error  in  the 
reconstruction  can  only  decrease  (or  stay  the  same)  at  each  iteration 
(14:197).  In  symbolic  form,  this  property  of  the  error-reduction 
algorithm  may  be  written  as 


<  EaK  < 


(146) 


2  2 

where  Ep  ^  and  E^  ^  are  the  errors  in  the  Fourier  modulus  and  object 
at  the  iteration,  respectively. 

From  the  form  of  the  equation  (146),  one  can  see  that  if  the  minimum 

error  in  the  Fourier  modulus  can  be  calculated,  and  if  at  the 

2  2  2 
iteration  Ep  ^  is  equal  to  the  minimum  error,  then  Ep  =  Ep  ^  . 

*  2  2  2 
Therefore,  it  would  follow  that  Eq  ^  =  %  k  w*ien  %  k  *s 

possible  value  of  the  Fourier  modulus  error.  In  this  manner,  then,  one 

arrives  at  the  least  amount  of  error  that  can  be  expected  in  the 

recons tructed  image. 

Calculation  of  the  Minimum  Error  in  Image  Reconstruction  by 
Phase  Retrieval  as  Applied  to  Transillumination. 

The  following  calculation  is  based  on  the  relationship  that  exists 


between  the  intensity  and  point  spread  function  (as  measured  in  the 
diffraction  plane)  to  the  Fourier  modulus  of  the  object.  This 
relationship  can  be  expressed  as 


(147) 


where  |r(f)j,  jl(f)|  ,  and  M(f)j  represent  the  magnitude  of  the  Fourier 
transforms  of  the  object  intensity,  image  intensity,  and  point  spread 
function,  respectively.  The  quantity  J M( f )  is  called  the  modulation 
transfer  function  (MTF).  The  frequency  f  is  measured  in  cycles  per  unit 
length. 

The  phase  retrieval  method  consists  of  making  an  estimate  of  the 
Fourier  modulus  of  the  object  by  taking  an  intensity  measurement  in  the 
diffraction  (Fourier  transform)  plane  and  dividing  this  quantity  by  a 
measurement  (or  estimate  based  on  a  theoretical  model)  of  the  MTF. 

This  estimate  of  the  Fourier  modulus  of  the  object  is  then  inserted 
into  the  error-reduction  algorithm  where  object  domain  and  Fourier 
domain  constraints  are  applied  to  arrive  at  a  new  estimate  of  the 
object.  The  steps  involved  in  the  algorithms  as  described  by  Fienup 
are 

"(1)  Fourier  transform  an  estimate  the  object;  (2)  replace  the 
modulus  of  the  resulting  computed  Fourier  transform  with  the 
measured  Fourier  modulus  to  form  an  estimate  of  the  Fourier 
transform;  (3)  inverse  Fourier  transform  the  estimate  of  the 
Fourier  transform;  and  (4)  replace  the  modulus  of  the  resulting 
computed  image  with  the  measured  modulus  to  form  a  new  estimate 
of  the  object."  (14:194) 

For  further  discussion  of  the  steps  of  the  algorithm  and  a  description 
of  the  function  (object)  and  Fourier  constraints  on  the  problem  the 
reader  is  referred  to  Fienup  (14:217). 

It  can  be  seen  that,  in  the  case  of  transillumination,  the  measured 
value  of  the  Fourier  modulus  of  the  object  would  be  corrupted  by  noise 
generated  by  the  scattering  process.  The  noise  in  the  Fourier  modulus 
F(f)|  ultimately  limits  the  efficacy  of  the  phase  retrieval  method 
(14).  Mathematically,  this  noise  is  expressed  in  the  point  spread 
function  or  its  Fourier  transform,  the  MTF.  In  practice,  the  MTF  could 
be  obtained  by  taking  a  measurement  on  a  reference  point  source  through 


the  scattering  medium.  Alternatively,  the  NfTF  could  be  derived  from  a 
calculated  point  spread  function.  (14:217). 

Since  the  point  spread  function  for  a  plane-parallel  mediun  has  already 
been  obtained  in  Chapter  II,  the  ITTF  can  now  be  calculated  and 
substituted  into  equation  (147)  to  find  the  value  of  j F( f )|  which  would 
be  used  in  phase  retrieval  as  the  measured  Fourier  modulus. 

Calculation  of  the  ffTF. 

The  reader  will  recall  from  Chapter  II  that  the  point  spread 
function  P(p)  may  be  written  as 

t?(p)  -  Pcco)<-  VP)  <148> 

where  P  <p)  is  the  coherent  (average)  intensity  and  P^(/9)is  the 
incoherent  (fluctuating)  intensity. 

From  equation  (148)  one  obtains 

IcfWcCO+FiCf) 

=  M(£) 

=  Mc(£)  +  MiC-H  (149) 

Equation  (149)  is  the  Fourier  transform  of  P^O)  and  M(f)  is  the  optical 
transfer  function.  The  quantity  M  (f)  represents  what  will  be  called  the 
coherent  portion  of  the  transform  and  M^(f)  the  incoherent  portion. 

At  the  outset,  the  non-normalized  incoherent  M^'(f)  and 
non-normalized  coherent  M  ' (f)  will  be  calculated.  Following  these 
calculations,  the  normalization  factor  will  be  obtained  to  arrive  at  the 
normalized  MTF. 

First,  M^(f)  is  the  Fourier  transform  of  the  incoherent  portion  of 

P  ph 

Mi  &)  = 


?i(£)ex^Ci27r£*p  )<3p 


(150) 


»Vl 


Pi  (p)  =  Io(V/$)x«t  [-r(i-w„)-  (p/fyyj 


and  2Hf  =  (k/f^)  P^  ( f q  =  focal  length  of  lens  in  Figure  29). 
Let  /5j  =  P3  (C0 S<p'£  +  5i^<pfy  ) , 

/5=p(cos<p><-  sin(py) 


then, 


p.p  =  PP>  (cos<pcos<p+  sin<priin<p)  =  pp^ Co^(<p -<p) 


dp  =  p<3pcl<p 


Therefore,  from  equations  (150)  and  (151): 


tAj(  )=l  l?i  (  )exp[i^p^cos(^V)J P  <3/®  ^ 


o  o 


e*p[i  Lpo  costs’-  <p)]  Pi<jo)  pdp 


o  I'O 


=  fznJo  C^PPd)Ei(p)pdp 

"'o 


-  2rno(a/pij1eKp[-T(l-K)] 

xA-(f0^)^f[-wjdp 


Now  using  the  formula 


f°° 

/  exp  [-  a?t2  (^ti)]  J0  Cbt)dt  = 


('Re*'>-1  ;  "Re  a?  >0  ) 


equation  (153)  becomes 


ifo )  r  7rl  p  a2  e  x  $  [- r  ( 1-  wc)]  exy  ^  ] 


(154) 


Now  to  calculate  M  '  (Pd)  (the  prime  indicated  denotes 


non-normalized).  First,  recall  that 


?ctp)  =  Ioexp(-r)Ca-/P')2Ji2(^a./f0; 


(155) 


and  let  b  s  ka/f^,  then 


W  -  I0expc-T)a?/  ^  Jt2 (b/5)e%p^(f(/f)p-pd ] dp 


(156) 


Following  steps  similar  to  those  that  led  to  equation  (152)  one  obtains 

r?“ 

Mc'(Pa)  =Iee"  &Jj  )ex£pVf/>^cos(<//-<p)pdpd<p 

Interchanging  orders  of  integration  results  in 

M'(  Pai  =  V*f(-  t)  a 2f 27ajK/WPd  (i/p)  J2  (bp)  dp 


Now  using  the  formula 


b  CU+v-X  +  1) 
z 


2Aa‘ 


y-\+t 


r^fDc 


N 

2  ' 


Q^(^+^-X+l)>  0;”Re\  >-l;Oib$a-] 


(157) 


V 


Setting  a=b,  fl  =  V  =  1,  Mc<Pd)  can  be  written  as 
M'(pa1  =  27Tl0exp(-T)a2  J0(Mpd) 


x  unto  . 


'Z^'- j.  (C  j  z) 


=  27Tl0ex|)  C-r )  a2  J0[(Vfoyy] 

v  i  rc&-c-<3) 

Src2>  rce-orce-d) 


=  27ric,exp(-r^a2Je[(^4JVdl 

x  i  r^Fco 

zr^  r^r<2> 


and,  therefore, 

M'cpj)  =  TTl^exf^-na2  Jo^VOPj)] 


(158) 


Finally,  combining  equations  (154)  and  (158)  and  noting  equation  (149): 

M'CA)  =  TDe^C-T)^le  (^=  A) 

^  0 

-  r\  r\  .  2.) 

4I0a  exppTd'W^Jextf 


(159) 


Equation  (159)  is  the  desired  result  for  the  non-normalized  WTF.  The 
first  term  on  the  right  side  of  equation  (159)  represents  the  coherent 
(average)  part  of  the  NTTF  and  the  second  term  represents  the  incoherent 
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(fluctuating)  part. 


To  get  the  normalized  MTF,  it  is  necessary  o  calculate  the  Fourier 
transform  of  the  point  spread  function  at  zero  frequency: 


Normal  i*z.  rng  Fax  “for  =  '(p )  dp 


=  j  \(P)+\(P)  <JP 


(160) 


2. 

=  I0C^')  exf  [-rd-w£,)](;l/2)Pi  Cn) 

1  ±  (161) 


Equation  (161)  is  the  incoherent  portion  of  the  normalization  factor. 
Now  set 


■  oo 


-  IcexpC-na2- .  I  jpl 


P 


(162) 


Using  equation  (157)  and  letting  c  =  ka/f^: 


•  oo 


?(p)^p~  I  e*fe(-T)az  [ — c  ^(1/2^  1 


4I>)  rc^/z) 

X  2J1  (l/z.;  -1/2.3  2;  1) 


-  Icexf  (-r)aL2c  Tl^2!.  — 

^  4T(3/2.)  TWW) 


-  I0ex£(-r)a*c(i/z)^ 

=  Ioex^(-r)azc.^_ 


Therefore, 


- 


S^dp-I^KVi3 ^  +  I0f^r  aze^[-r(l-W0)] 


Now  dividing  M  1  (Pd)  +  Nh'(Pd)  by  the  normalizing  factor  one  gets  for 


the  optical  transfer  function  /1C/D  D.  Y" 

M (pd)  -7raz  J0(4-Pl)  +-  7razexf  (rW^ex^  -  - 1  -1  - 


(16' 

and  the  MTF  =  M<Pd)  . 

From  equation  (147): 

N  ■  ^ 

|l(/¥>l 

Equation  (165)  may  be  rewritten  as 

%  = 

IX<P4'i 

1+  |MiQ) 

(16i 

It  is  clear  that  Jl(Pd) j/ |Mc(Pd)|  is  the  correct  Fourier  transform  of 

the  object.  Denoting  by 

F(pd)lthe  exact  modulus  of  Fourier  transform 

of  the  object,  equation  (166)  becomes 

d  1+ 

|WC^ 


(167) 


Subtracting  |F(pd)|  from  both  sides  of  equation  (167)  and  squaring  and 
dividing  the  result  by  |F(Pd)  ^  yields 

l^l2  1  4  IWI 


|IAc<Pa)l 

where  the  left  side  of  the  equation  is  the  definition  which  the 
error-reduction  algorithm  gives  the  error  in  the  Fourier  modulus 
estimate.  Therefore, 


1+  I 

|Mc«V>l 

Now  if  one  can  find  the  minimum  value  of  Nb(Pd)  /  (Pd )  ,  then 
substituting  that  minimum  value  into  equation  (168)  will  yield  the 
minimum  error  that  can  be  expected  in  the  estimate  of  the  Fourier 
modulus. 

The  ratio  NT  (Pd)  /  McPd)|  is 

r  [«/aW 

[  Nlj(  j)  _  e^(r Wo)e^tr  4 


(168) 


(169) 


This  equation  has  many  relative  minima.  In  fact,  the  absolute  minimum  of 
equation  (169)  occurs  at  Pd  =00.  At  Pd  =00,  equation  (169)  implies  that 
the  error  in  the  estimate  of  the  Fourier  modulus  is  zero  (by  equation 
(168)).  However,  there  is  a  practical  limnit  on  the  magnitude  of  Pd 
since  Pd  has  a  cutoff  value  of  2a  (the  diameter  of  the  aperture). 

Without  further  argumentation,  therefore,  it  is  difficult  to  assess  the 
implications  of  equation  (169). 

Consider,  however,  the  following  argument  which  will  show  that,  for 
very  large  T,  the  only  way  to  obtain  a  mathematically  defined  estimate 
of  the  Fourier  modulus  is  to  evaluate  Pd  at  approximately  its  zero 
value. 

First,  evaluate  M^Pd)  +  M^(pd): 


n£j0  (jPj)  +ta.lex-f  (TW0)ex-£ 


-(£/W 

4 


.Hi  + 

fo37T 


(7T) 

ZPl 


-a.2  ex-^CrWP) 


For  largeT(and  any  value  ofPd),  this  expression  becomes 

-(kp.P-'f 

......  . .  .  ,  7 r  a2  ex t»(7~Wg>x£  -  /  — 

ey^CrWo) 


(i7o; 


Uc(fi^^i(pd)  -  J 


Recalling  that 


nz  -  c  K.  i-  rTWo't 

p-i  -  (**10 


equation  (170)  becomes  ,,  ,  ,  -j  rZ  \k,  XT  *~ 

L  lCnr  i  tvMV 

Mito  -  2mUZ  3^JS 


=  2(TD1/Piext,(  ^WeKr 
1  M  4«c 


Now  using  the  same  values  forT,  Wq,  kf  andOk'  that  were  used  in  Chapte 


It  can  be  seen  that  for  any  value  of  P^  (0  <  P^  <  2a)  other  than  zero, 
equation  (173)  is  virtually  zero  due  to  the  exponential  factor. 


Consequently,  by  equations  (170)  and  (171)  the  Fourier  modulus  estimate 


becomes  undefined.  Therefore,  to  get  a  meaningful  Fourier  modulus 


r»fiv#vr.  '* 


estimate,  /O^  must  be  evaluated  at  virtually  its  zero  value. 
Substituting  P^  =  0  into  equation  (164)  one  gets 

|‘Mi(o)|y|Mc(rt|  =  ex-pCrW.) 

Now  substituting  equation  (168)  the  result  is 

E  2  =  I  1/(1+  ex-p^Vll2 


(174) 


(175) 


For  large  T  and  Wq~  1,  Ep  =  1,  or  in  other  words,  the  error  in  Fourier 
modulus  estimate  is  virtually  100  percent. 

Equation  (175)  expresses  the  minimum  error  one  can  expect  in  the 
estimate  of  the  Fourier  modulus.  Therefore,  by  the  unequality  relation 

EF,ku  *  ^O/k 


characteristic  of  the  error-reduction  algorithm, 

E02  =  E l  =  1 1/C  l  +  ex-p  (T  W0)>  1  |x 


(176) 


for  large 7".  Consequently,  it  can  be  concluded  that  the  error  in  the 
reconstructed  image  is,  for  all  intents  and  purposes,  100  percent  and 
that  the  phase  of  the  object  cannot  be  retrieved. 
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Appendix  I:  The  Scattering  Cross  Section 


The  scattering  cross  section  (j  may  be  defined  in  terms  of  the 

s 

differential  scattering  cross  section  (J^  where 

(J-3  (6,  l )  =  lim  [(^  SsVSi] 

d  K-co 

=  (177 ) 

R  is  the  distance  from  the  scatterer  to  the  point  of  observation,  and 
Sg  are  the  magnitude  of  the  incident  and  scattering  power  flux  density 
vectors,  a  is  the  total  scattering  cross  section  (O’  +  O" )  and  pC^,*) 

L  S3 

is  the  phase  function. 

The  total  observed  scattered  power  at  all  angles  is  produced  by  a 
particle  with  scattering  cross  section  CTS' 

°s=  [%<& 

->4 n 

=  [ 

JtfT 

(178) 

where  dCJis  the  differential  solid  angle. 

The  albedo  Wn  of  a  particle  is 


(178) 


w„  = 


-  /  V(6, £)<3W 
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(179) 


vjNJoW .vA‘0. 


Appendix  J:  Solution  to  the  Transport  Equation  Using,  the 


Large  Particle  Approximation  (19:236 


As  shown  in  Chapter  II,  if  the  size  of  scatterers  is  much  greater 
than  a  wavelength,  it  is  possible  to  simplify  the  equation  of  transfer 
(transport  equation)  using  the  small  angle  approximation  (see  Chapter 
II,  page  ).  One  then  obtains 

s-Vti^p,S) 

O' z. 


=  f -pcs-  ?')!(*,  P,  ' 

47T  JJ 


(180) 


where  r  =  xx  +  yy  +  zz  =  O  +  zz,  +  %-y, 

'  £■  /“it  / 


£  =  lx  +  my  (19:236).  y 

A  general  solution  of  equation  (180)  may  be  obtained  by  a  Fourier 
transform  method.  Denote  by  1^  the  Fourier  transform  of  I  with  respect 
tOy0(19:236): 

TifeAs)-  /lftAs)ex  (i8i) 


W,s)  =  lW)exf(-i*-P)tf 


(182) 


Now  the  Fourier  transform  of  equation  (180)  may  be  obtained.  First, 


note  that 


=feW,s) 


(183) 


which  is  the  Fourier  transform  of  the  first  term  on  the  left  side  of 


equation  (180)  (19:236). 

To  get  the  Fourier  transform  of  the  second  term  note  that 

s-VtI(?  .?)>  jSM+ffa ifc,  ,i)y] 


(184) 


K'p  -  (Kxx+K,9).(*Si*y$) 


-  *yX  +-  K y 

*  J 


and,  therefore, 

i-V,  iM,?)  = 


-it*  ~S2.  /  ^'^/c 

<z7r>4  * 

.  \  f  T  A  -  i  AC  ,/5 

-1S*  imtW  1  y^e 

1  J-co  r°°  .  x 


Jroo  .  -»  ^J~00 

rr1e-1<,pcl/c 

-od 

r°°  _iK.p 

r  (OT)M  ^(185) 

-oo 


-lAC-p. 
Ae  cU 


:$a 

iyq 


But  equation  (184)  expresses  the  Fourier  transform  relation: 


(186) 


which  yields,  then,  the  Fourier  transform  of  the  second  term  on  the  left 
side  of  equation  (180). 

The  Fourier  transform  of  the  first  term  on  the  right  side  of 
equation  (B-l)  is  straightforward: 


(187) 


The  Fourier  transform  of  the  second  term  on  the  right  side  of 
equation  (180)  is  oo 


s(^,f , S)ds'  • 


r  Wt 


Oo 


±7 T 


JJvCS-s'rfflCifi,  s)}ds' 


-  OO 


oo 


=  ^  f  ffcs-f) Let/,  S)di'  <188> 

4tt  J  I*  x 


-oo 


Collecting  equations  (183),  (186),  (187),  and  (188)  one  obtains 


oo 


(189) 


-oo  r  0 

It  is  convenient  to  simplify  the  first  term  of  equation  (189)  by 


letting  (19:237) 

ItCz., s)=  IzCe/  *•] 


(190) 


then  equation  (189)  becomes 

oo 


^-^[JpcS-S'Jexfjd  -cS-S'jzjl^A,^' 


=  O 


Now  using  the  relations^ 

=[h(i)exf(iS‘i)iS 

(192 

J  J-OO  r  f-oo 

=1  k2  (Z,*,  S>x£(i5-‘i>i  t  (193 

J  J-Oo 

the  Fourier  transform  of  the  convolution  integral  in  equation  (191)  is 
taken.  First,  let  r  =  s  -  s'  and  dr  =  ds  then, 

-  f  r^x^ex-gC-i^-rjeAf 


=.m 


(194; 


=  -  mj 

-Oo 

Substituting  equation  (194)  into  equation  (191)  yields 

[ff  [[ f(r)eir'(*~*^dr\l  *, ds ’ 

JJ[JJ  > 


=  pc«r-<fe)F(*,*.q> 


(by  equations  (192)  and  (193)). 


The  Fourier  transform  of  the  first  term  in  equation  (191)  is 
straightforward:  ^ 

A  h:1^'  = S)e  ^ t 


.vVfVv«! 


«V  ‘  ‘Avi*. 


v,;*\  y 


y^wwvSSSScs 


*Vi 


Then,  one  obtains  from  equations  (195)  and  (196): 


(197) 


Equation  (197)  is  exact  and  may  be  solved  using  the  integration  factor 

eKf  J  ft&k  : 

exp  jj BP<k‘  •  ^  +-"BPexp  (19g 


where 


-b=/*Ql  ,?  =  K3-*0  r  =  Ffe,<|). 

Ar7T  / 


Therefore, 


4_  [exp(fvPd£)Il  =  0 


e>{>[ /BP  Az'jFrC- 


(199) 


where  C  is  a  constant  of  integration.  Setting  z  =  0  the  constant  c  may 


be  evaluated: 


L  Jo  J  e-o 


C=  FCo,x,  5) 


From  equation  (199)  one  then  obtains 


expQ  /  w&'3FM,?>=:F(o,M) 


***N>*!«y*g 


1  r<  rl 


(200) 


The  final  solution  I(z,p,'s)  for  the  small  angle  approximation  may  now  be 
derived.  First,  ^ 

1M,s)=(4  /  Il<'Z/.s)exfflK-p)d/f 

-Oo 

xex-f  (201) 


smce  1  (■a,xJS>=  -i-  /  F(z,<  1)e'1S^<l4 

<OT>  4 

equation  (201)  becomes  qq 

Hz  A  3)  -  JL*  y"dxexp^-i?.^+ 

X  faeifC-iS'UFfr/A) 


(202) 


Now  using  equation  (200)  and  changing  the  variable  q *  =  q  -Kz  and  then 
removing  the  prime  (q  ->  q  +AT z): 

I  (Z,p, s)  =  i_4^Scexp|i<-/?+ 


x  Jd<exf(-iS-?>F0(  A^?[jjr-P(q-  ] 


-  I 


C27r)^  / 

V  -i-00  /-OO 

X  /clqexf[-iS.(q4.t*3F6(x^+?c«) 

JLco  rZ 


rM’»  r*.  •*« 


Then,  simplifying  terms  I(z,p,s)  becomes  _oO 


xe*fH /Vtt1-  |^S+*)<z (203) 

Equation  (203)  is  Ishimaru's  solution  for  I(z,p,s)  (19:238). 
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Appendix  K:  Calculation  of  the  Fourier  Transform  of  the 

Gaussian  Phase  Function 


The  phase  function  describing  the  angular  distribution  of  scattered 


radiation  can  be  written  as 

£(  5}  =  4.0^  Wo 


(204) 


To  find  the  Fourier  transfoii:  P(«J)  of  p(S)  the  following  formula  is 


used: 


O 0 


'  4a^W e>exf(-°{£  sx)exf(is.  £ )  ds 


(205) 


Letting 


-DO 


5=  sin#(cosv?x 

q  -  (coyfa  +-  sinewy  ) 

then, 

5.^  -  <qsiT)0cos^cos(^-f- 

~  SXTl  6  |jC05^COS(^,f'+-  SXTl  ^Sl'H 

=  c^sin  9  CcOS 

(206) 


Substituting  into  equation  (205): 

.00 


V(^)  -  I  [4a^V/0exjb 8 [cos  (y>  -  <//)]  j  I5 


-CO 


(207) 


Rearranging,  and  changing  order  of  integration  gives 

P (q_)  =  J 4 .Qfe % e*1? C-Q£ )  j  J  exf> Oco s 

X  cos  8  sin#  d# 

rn 

~  J  4<XpW027rJ0(qsiri^)exf  o^,  sin2# ) 

X  cos#siri$<A# 

rn 

=  4ttW0  /  ^  6jsin0)ex$(-Q^sinz0)2a^  ( 

Since  .  /j°  x cos# sin#  d# 

Since  IS  I  =  sxni95 

P<q}  =  47rW0/j0(q.^ex^C-<^Sl)2A^scls 
•'o 

=  frTWo  (zatf)  I  J0(qs^ex^C-a^s2)5cls 

roo 

=  47rw0(2a^jf  ^(q_s)ex^[-[(^y^2slJsis  < 

Now  use  the  formula  (1) 

rOo 

J  J^Cat)  c  x  £  (-p2  t2-)  1  At 

to  get  ~  ^J?5meXf£f£‘) 

P(^  =  4r7TW0  czotf).  X- •  e*f  (- 

-  47rWfl  exp  (-  aX/4 Ap) 


to  get 


which  is  the  desired  result. 
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